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USERS MANUAL FOR THE VARIABLE DIMENSION AUTOMATIC SYNTHESIS 


PROGRAM (VASP) 

John S. White and Homer Q. Lee 
Ames Research Center 


SUMMARY 


VASP is a Kariable dimension Fortran version of the Automatic Synthesis Program, a computer 
implementation of the Kalman filtering and control theory. It consists of 3 1 subprograms for analy- 
sis, synthesis, and optimization of complicated high-order time-variant problems associated with 
modern control theory. These subprograms include operations of matrix algebra, computation of 
the exponential of a matrix and its convolution integral, solution of the matrix Ricatti equation, and 
computation of dynamical response of a high-order system. 

Since VASP is programmed in Fortran, the user has at his disposal not only the VASP 
subprograms, but all Fortran built-in functions and any other programs written in the Fortran 
language. All the storage allocation is controlled by the user so the largest system that the program 
can handle is limited only by the size of the computer, the complexity of the problem, and the 
ingenuity of the user. No accuracy was lost in converting the original machine language program to 
Fortran. 

The principal part of this report contains a VASP dictionary and some example problems. The 
dictionary contains a description of each subroutine and instructions on its use. The example prob- 
lems give the user a better perspective on the use of VASP for solving problems in modern control 
theory. These example problems include dynamical response, optimal control gain, solution of the 
sampled data matrix Ricatti equation, matrix decomposition, and pseudo inverse of a matrix. 

Listings of all subroutines are also included. 

The VASP program has been further adapted to run in the conversational mode on the Ames 
360/67 computer. The necessary procedures are given in appendix C. 


INTRODUCTION 


The VASP, the Variable dimension Fortran version of the Automatic Synthesis .Program, is the 
new Fortran IV version of the ASP, the Automatic Synthesis Program. It is intended to implement 
the Kalman filtering and control theory. Basically, it consists of 3 1 subprograms for solving most 
modern control problems in linear, time-variant (or time-invariant) control systems. These subpro- 
grams include operations of matrix algebra, computation of the exponential of a matrix and its con- 
volution integral, and the solution of the matrix Riccati equation. The user calls these subprograms 
by means of a FORTRAN main program, and so can easily obtain solutions to most general prob- 
lems of extremization of a quadratic functional of the state of the linear dynamical system. 



Particularly, these problems include the synthesis of the Kalman filter gains and of the optimal 
feedback gains for minimization of a quadratic performance index. 

The VASP is an outgrowth of ASP, which was developed for NASA under contract with the 
Research Institute for Advanced Studies, a division of the Martin Company. There are two urgent 
reasons for reprogramming ASP into the present Fortran version. First, ASP was programmed in 
FAP (Fortran Assembly Program) and could be used only on the IBM 7090-7094. Second, many 
complicated time-variant analysis, synthesis, and optimization problems tax the capability of the 
ASP and other programs written in the Fortran language. Fortran IV language makes the program 
adaptable to a much wider class of computers and expands its versatility. 

I 

The VASP is based extensively on a Fortran version of ASP, called FASP (Fortran ASP) by its 
programmer Mr. Don Kesler of Northrop, Norair. 

Two basic questions the user will inevitably ask are: 

(1) How accurate is VASP compared with ASP? 

(2) What is the highest order of system that VASP can handle? 

The answer to these questions depends on the number of significant digits carried by the user’s 
computer and the amount of available storage in the computer. To answer the first question in a 
more concrete way, the check cases given in the ASP manual were duplicated and the results were 
compared with those in the manual. The accuracy of VASP was found to be the same as that of 
ASP. The second question can best be answered by first noting some of the basic differences 
between FASP and VASP. The pertinent difference between the two is that VASP has variable 
dimensioning and more efficient storage. To allow a computer to handle the highest order system 
possible, all matrix storage is assigned by the user’s main program. Consequently, as an illustrative 
example, a 125,000-byte version of the IBM 360/50 can easily determine the solution of the matrix 
Riccati equation for a 30-order system (perhaps 40, depending on the size of other related matrices). 
Another basic difference between these two Fortran versions is that VASP diagnostics are more 
self-explanatory. 

To recapitulate, the objectives of VASP are flexibility and versatility so that it can serve the 
maximum number of users. To achieve these goals FASP was revised extensively so as to have, for 
example, variable dimensioning, more efficient storage, and more self-explanatory diagnostics. 

In this report, no attempt will be made to discuss any details of the theoretical background and 
the algorithms associated with the appropriate subprograms since they are well documented in the 
ASP manual, an NASA contractor report (NASA CR— 475, 1 966). This report does not repeat 
information from the contractor report, and the user is urged to consult that manual so that he may 
utilize VASP proficiently. 

This program can be obtained from the centralized facility known as COSMIC, located at the 
Computer Software Management and Information Center, Barrow Hall, University of Georgia, 
Athens, Georgia, 30601. 
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FEATURES OF THE PROGRAM 


The advantages of VASP over ASP are (1) a more versatile programming language, (2) a more 
convenient input/output format, (3) some new programs, and (4) variable dimensioning. 

Since VASP is programmed in Fortran, it can be used in a very large class of machines. 
Moreover, as VASP is part of a larger main program, all the Fortran built-in functions are available 
to the main program. Furthermore, any subroutine available in the Fortran language may be used. 
In other words, the user has at his disposal the combined capabilities of VASP, Fortran built-in 
functions, and all other subroutines written in Fortran. 

The input/output subroutines have been changed and now consist of READ, RDTITL, and 
PRNT. In addition, LNCNT has been added to control paging. The VASP allows the user the 
option of using the existing standard VASP format, or of supplying the output format of his own 
choice. For a more detailed explanation of how to exercise this latter option, see the dictionary 
entry under PRNT (p. 10), or Example 2. 

Our experience with ASP is that certain groups of statements are often repeated. For the 
user’s convenience, these groups of statements are incorporated as new subroutines in the VASP. 
They are AUG, UNITY, and SCALE. Detailed explanations of them are available in the VASP 
dictionary in this report. 

To utilize the storage space as efficiently as possible, the subroutines are written with variable 
dimensioning, with the storage allocation controlled by the user’s calling program. Consequently, 
it is necessary to provide some dummy storage space as a part of the argument of the subroutine. 
From the user’s point of view, the price for efficient storage is inconvenience. All the subroutines 
are written in double precision for adequate accuracy; that is, all matrix and scalar variables , except 
integers, are assumed to be in double precision. 


Universal Features 

The arguments in the subroutines are arranged in the following order: 

Input arguments 
Output arguments 
Dummy arguments 

These are also arranged so that matrices occur before scalars. 

An array of length two must be allocated by the user to store the dimensions of the matrix, 
and this array must be included in the subroutine call statements. For example, the add subroutine 
is called by 

| . 

CALL ADD(A,NA,B,NB,C,NC) 


and performs the matrix operation 
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C = A + B 


Here NA, NB, and NC are arrays of length two which contain the dimensions of matrices A, B, C, 
respectively. In other words, the numbers of rows and columns of matrices A, B, and C are stored in 
NA, NB, and NC, respectively. Specifically, the number of rows of A is stored in NA(l)andthe 
number of columns of A, in NA(2). 

In general, the dimension array associated with an input matrix contains input information to 
the subroutine, while that associated with an output matrix contains output information. The dic- 
tionary shows the few cases where this rule is violated. In the example above, dimension arrays NA 
and NB 1 are inputs (since matrices A, B are inputs) and must be loaded before entering this sub- 
routine. On the other hand, NC is an output, since C is an output. That is, the values of NC(1) 
and NC(2) are computed in the subroutine and are available to the calling program upon return. 

When a dummy array is required, it must be appropriately dimensioned in the calling program. 
The required size is given in the appropriate dictionary entry. 

Most of the routines check the array dimensions for compatibility and reasonableness, and 
check for adequate storage available in the DUMMY array. The “reasonableness” test is to see 
that all matrix dimensions are greater than zero, and that the product of the matrix dimensions is 
less than the constant MAXRC. In the program MAXRC has been set at 6400. It is recommended, 
however, that the user reset MAXRC to equal the size of his matrices, and thus prevent accidental 
overflowing of the matrix dimensions. If the matrices are incompatible or unreasonable, or if a 
mathematical error is found, a self-explanatory error message is printed, and no further computations 
are made in that subroutine. However, computation does go on to the subsequent steps, which are 
likely to be wrong also. After 1 0 such errors, the program is terminated. 

The VASP program uses double-precision arrays, so that the user’s main program must define 
each matrix to be double precision, and to contain a sufficient number of cells to accommodate the 
matrix. The dimension statement may classify the array as one- or two-dimensional, or even more. 

For example, to use the matrix A, which is a 5 X 5 matrix, any of the following dimension 
statements will be adequate: 

DOUBLE PRECISION A(25) 

DOUBLE PRECISION A(5,5) 

DOUBLE PRECISION A(3,10) 

DOUBLE PRECISION A(100) 

The important factor is the total number of cells reserved, and the user may reserve more cells than 
the matrix requires, or, conversely, he may put a smaller matrix than originally planned in a specific 
array. The VASP program stores data in an array as a string of columns, just as Fortran does. 


1 The convention used here, and throughout this report, is that the name of a dimension array is obtained by 
prefixing the letter N to the matrix name. 
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However, it stores the matrix A according to the dimension given in NA, whereas Fortran stores 
A according to the dimensions in the Fortran dimension statement. 2 

Consider the following example. The Fortran statements are: 

DOUBLE PRECISION A(5,5), B(5,5), C(5,5) 

DIMENSION NA(2), NB(2), NC(2) 

CALL READ(3,A,NA,B,NB,C,NC, . . . ) 

The first card in the data deck specifies NA = 5,5, followed by cards with 25 data words for A; then 
NB = 4,4, followed by 16 data words of B; finally, NC = 6,6, followed by 36 data words of C. Since 
the storage of data in VASP is controlled by the VASP dimensioning, the 25 words for A will 
exactly fill the reserved storage and the 1 6 words for B will fill the first 1 6 cells of the storage 
reserved for B. The 36 words for C will completely fill the reserved storage for C and overflow 
into something else. The user can prevent this overflow by setting the test constant MAXRC to 25. 
The error test in the READ subroutine will note that the product of NC( 1 ) and NC(2) is greater than 
MAXRC, and will return an error message. This selection of MAXRC will limit all other VASP arrays 
to 25, so it is frequently desirable to dimension all arrays the same. 

Occasionally the user may wish to refer to a single element of a matrix. Since FORTRAN 
statements use the FORTRAN dimension statement, a reference to B(4,4) in the previous example 
will select the 19 element in the B storage. However, VASP, using the VASP dimension, has stored 
B(4,4) in the 1 6 element of B, which is not the same. To actually select a specific element, say 
B(i,j), one must refer to B((j— l)*NB(2)+i,l). In the above example, references to A(ij) will be 
correct, since the FORTRAN and VASP dimensions are the same. 


System-Dependent Features 

Two subroutines in the VASP package are system dependent. The first is BLKDTA. Data 
statements in this subroutine control the printing. They require a printer with at least 1 1 5 charac- 
ters per line, and place 45 lines on each page. These requirements may be changed as needed. The 
second is ASPERR, which calls a system subroutine for error tracing. The description of ASPERR 
indicates any necessary changes to match the system. 

The VASP programs frequently generate very small numbers. The computer operating system 
may detect these small numbers as underflows, and print error messages. If so, the user should 
arrange to turn off the underflow error messages. 


THE VASP DICTIONARY 


A detailed description of all the subroutines is given in this dictionary. Each entry is 
organized into subheadings that describe the subroutine and explain how to use it. Other 


2 The storage in VASP is also compatible with the storage of “general” matrices in the IBM scientific subroutine 
package. 
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subheadings, such as motivation and remarks, are sometimes added to offer the user a better 
understanding of the theoretical background of the subroutine. 3 

The dictionary proper lists only those routines that the user is expected to call directly. Several 
additional subroutines, used internally, are also a part of VASP. The user may, however, wish to call 
these routines himself, since they are quite flexible. These additional routines are described in 
appendix A, and a complete listing of all programs is given in appendix B. 

Several procedures have been written to facilitate the use of VASP on Ames time-share system. 
Their usage and listings are given in appendix C. 

Table I lists all subroutines with their calling sequence, and the TSS procedures, for easy 
reference, while table 2 lists the approximate storage used by each of the VASP subroutines when 
compiled on the NASA Ames 360/67, OS system. Table 2 also lists the external references for each 
subroutine. 


3 Some of the subroutines are almost direct copies from the Northrop FASP. The detailed description of the 
theory is either obvious, or is given in the ASP manual (NASA CR— 475). Other routines were written by one of the 
authors. These were quite simple, and needed little description. Subroutines ANDRA, BDNRM, DECOM, and PSEU 
were written by John Andrews, Information Systems Company. Since no description of these subroutines is available 
elsewhere, a detailed description of their theory and usage is included. Because there were various programmers, the 
nomenclature internal to the various subroutines is not completely consistent. 
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TABLE 1.- SUBROUTINE CALL STATEMENTS IN VASP 


1. CALL READ(I,A,NA,B,NB,C,NC,D,ND,E,NE) 

2. CALL RDTITL 

3. CALL PRNT(A,NA,NAM,1P) 

4. CALL LNCNT(N) 

5. CALL ADD(A ) NA,B,NB,C,NC) 

6. CALL SUBT(A,NA,B,NB,C,NC) 

7. CALL MULT(A,NA,B,NB,C,NC) 

8. CALL SCALE(A,NA,B,NB,S) 

9. CALL TRANP(A,NA,B,NB) 

10. CALL INV(A,NA,DET,DUM) 

11. CALL NORM(A,NA,AN0RM) 

12. CALL UNITY(A,NA) 

13. CALL TRCE(A,NA,TR) 

14. CALL EQUATE(A,NA,B,NB) 

15. CALL JUXTC(A,NA,B,NB,C,NC) 

16. CALL JUXTR(A,NA,B,NB,C,NC) 

17. CALL EAT(A,NA,TT,B,NB,C,NC,DUMMY,KDUM) 

18. CALL ETPHI(A,NA,TT,B,NB,DUMMY,KX)UM) 

19. CALL AUG(F,NF,G,NG,RI,NRI,H,NH,Q,NQ,C,NC,Z,NZ,II) 

20. CALL RICAT(PHI,NPHI,C,NC,NCONT,K,NK,PT,NPT,DUM,KDUM) 

21. CALL SAMPL(PHI,NPHI,H,NH,Q,NQ,R,NR,P,NP,K,NK,NCONT,DUM,KDUM) 

22. CALL TRNS1(F,NF,G ) NG,J,NJ > R,NR ) K > NK,H,NH,X,NX,T,DUMMY,KDUM) 

23. CALL PSEUDO(A,NA,B,NB,DUM,KDUM) 

24. CALL DECGEN(R,NR,G,NG,H,NH,DUM,KDUM) 

24a. CALL DECSYM(R,NR,G,NG ) H,NH,DUM,KDUM) 

Programs 25 through 31 are called internally and need not be used by the programmer. They 

are described in appendix A. 

25. CALL READ1 (A,N A,NZ,N AM) 

26. CALL ASPERR 

27. BLKDATA (nonexecutable) 

28. CALL PSEU(A,B)C,EE,DEP,IP,D) 

28a. CALL PSEUP(A,B,C,EE,DEP,IP,D) 

29. FUNCTION BDNRM(NR,CT,EE,D,KRV) 

29a. CALL TTRM(NR,CT,EE) 

30. CALL ANDRA(B,T,DPR,JP) 

31. CALL DECOM(A,B,C,E,JL,DCM,KP,D) 

The remainder of the items are procedures to facilitate the use of VASP on the Ames TSS. 

32. VASP$$ [inputdsname] [.outputdsname] 

33. CHNGIN [inputdsname] 

34. CHNGOUT [outputdsname] 

35. CMPL 

36. CLRVASP 

37. CONVASP [matrixsize] [,$A=name] [,$B=name] [,$C=name] [,$W=name] [,$X=name] [,$Y=name] [,SZ=name] 

38. RECMPT 

39. REWRT [n] 
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TABLE 2.- APPROXIMATE STORAGE REQUIREMENTS AND EXTERNAL REFERENCES 


Subroutines 

Storage 
decimal bytes 

External references 

1. READ 

1000 

READl , PRNT* 

2. RDTITL 

400 

LNCNT 

3. PRNT 

1400 

* 

4. LNCNT 

500 

None 

5. ADD 

800 

* 

6. SUBT 

800 

* 

7. MULT 

1100 

* 

8. SCALE 

700 

* 

9. TRANP 

800 

* 

10. INV 

2500 

* 

11. NORM 

1000 

* 

12. UNITY 

700 

SCALE* 

13. TRACE 

700 

* 

14. EQUATE 

700 

* 

15. JUXTC 

1000 

* 

16. JUXTR 

1100 

* 

17. EAT 

3200 

ADD, MULT, SCALE, NORM, UNITY, EQUATE* 

18. ETPHI 

2300 

ADD, MULT, SCALE, NORM, UNITY, EQUATE* 

19. AUG 

3300 

MULT, TRANP, EQUATE* 

20. RICAT 

5100 

ADD, MULT, INV, EQUATE, PRNT* 

21. SAMPL 

3700 

ADD, SUBT, MULT, TRANP, PSEUDO, PSEU, BDNRM, 



ANDRA, PRNT* 

22. TRNSI 

5000 

ADD, EAT, SUBT, MULT, PRNT, EQUATE* 

23. PSEUDO 

900 

PSEU, BDNRM, ANDRA* 

24. DECGEN 

2600 

MULT, TRANP, INV, NORM, EQUATE, DECOM, PSEU, 



BDNRM, ANDRA* 

25. RE ADI 

800 

PRNT* 

26. ASPERR 

400 

None 

27. BLKDATA 

None 

None 

28. PSEU 

5900 

MULT, NORM, BDNRM, ANDRA* 

29. BDNRM 

1500 

MULT, NORM 

30. ANDRA 

2000 

LNCNT 

31. DECOM . 

1500 

MULT, NORM, PSEU, BDNRM, ANDRA* 

COMMON/MAX/ I 



COMMON/LINES/ } 

200 


COMMON/FORM/ j 



TOTAL 

53,600 



*LNCNT and ASPERR are additional external references. 





1 . 


READ 


DESCRIPTION 

This subroutine reads 1 to 5 matrices from cards, along with the names and dimensions, and 
prints the same information. For each matrix the routine first reads a header card containing a four- 
character title, followed by two integers giving the row and column size of the matrix, using format 
(A 4, 4x, 214). Then the matrix data are read using READ1, each row of the matrix starting on a 
new card, using format (8F 10.2). If the card data is in exponential form, it must use a D exponent. 
The matrix title and the matrix are then printed using subroutine PRNT. 

If the header card contains no row and column size (i.e. , n = 0), then the matrix in storage is 
unchanged, no data cards are read for that matrix, and the previously stored matrix is printed. 

USAGE 

CALL READ(I,A,NA,B,NB,C,NC,D,ND,E,NE) 

Input Arguments 

Control constant: I 

where I is an integer from 1 to 5 and indicates the number of matrices to be read. If I is less than 
5, the extra matrices in the call list may be dummy variables, or repeated references to the same 
matrix; for example, 

CALL READ( 1 ,A,NA,A,NA, A,NA,A,NA, A,NA) 

Output Arguments 

Matrices: The first I of the matrices A,B,C,D,E 

Dimension arrays: The first I of the arrays NA,NB,NC,ND,NE 


2. RDTITL 
DESCRIPTION 

This subroutine reads a single card in hollerith format, and loads it into the array TITLES. It 
then calls LNCNT(IOO). This latter program in turn skips the printer to the next page, prints the 
hollerith information in the array TITLES, and positions the output to print next on line 3. 

USAGE 

CALL RDTITL 
It has no arguments. 
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3. 


PRNT 


DESCRIPTION 

This subroutine prints a single matrix, with or without a title line, and either on the same page 
or a new page. The matrix is printed using format (1P7D16.7) for the first line, and (3x,lP7D16.7) 
for all subsequent lines. The “3x” indents continuation lines for easier reading. 

REMARKS 

The standard format is stored in arrays FMT1 (for the first line) and FMT2 (for subsequent 
lines) both of which are stored in labeled COMMON as follows: 

COMMON/FORM/NEPR, FMT1 (6), FMT2(6) 

where NEPR is the number of columns of data to be printed (7, in the standard case). The standard 
format is loaded into COMMON in the BLKDATA program. If other formats are desired, they can 
be obtained either by changing the BLKDATA program, or by having the users main program change 
the contents of COMMON. 

CAUTION 

In writing a data statement for the formats, put FMT1 and FMT2 in separate statements, as in 
the BLKDATA program. If they are loaded in one statement, they will probably load incorrectly, 
because of the dimensionality of FMT1 and FMT2. Also NEPR must be consistent with the numbers 
in FMT1 and FMT2. 

USAGE 

CALL PRNT(A,NA,NAM,IP) 

Input Arguments 
Matrix: 

Dimension array: 

Matrix name: 


Control constant: 

2 heading, new page 

3 no heading, same page 

4 no heading, new page 


A 

NA 

NAM 

If NAM = 0, a blank name will be printed. NAM should contain four 
hollerith characters. It can be written in the calling sequence as 
4HAbbb. If written 1HA, the last three characters of the printed 
name will be garbage. 

IP= 1 heading, same page 


Output Arguments 
None 
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4. LNCNT 
DESCRIPTION 

This subroutine keeps track of the number of lines printed, and automatically pages the output 
as required. It is completely internal, and the user need not refer to it unless he has WRITE state- 
ments of his own. In that case, the user may (should) put the statement CALL LNCNT(N) before 
each WRITE statement, where N is the number of lines to be printed. 

Page length is controlled by the variable NLP set in the BLOCK DATA program to 45. This is 
an installation-dependent variable, and may be changed as necessary. 

The subroutine provides one line of print at the top of each page. This line contains 92 
characters, of which the first 72 are available for the programmers use and may be loaded by use of 
RDTITL. The remainder contain “VASP PROGRAM.” The 92 characters are contained in the 
array TITLES, which is, in turn, contained in the COMMON area LINES. If N > NLP, the printer 
will automatically skip to the top of the next page, and print the title line. 

USAGE 

CALL LNCNT(N) 

Input Arguments 

Constant N = number of lines to be printed 
Output Arguments 

None 


5. ADD 

DESCRIPTION 

This subroutine computes the matrix sum 

C = A + B 

USAGE 

CALL ADD(A,NA,B,NB,C,NC) 

Input Arguments 

Matrices: A,B 

Dimension arrays: NA,NB 

11 
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Output Arguments 


Matrices: C 

Dimension array: NC 

REMARK 

Matrices A and C may share the same storage space or matrices B and C may share the same 
storage space. 


6. SUBT 
DESCRIPTION 

This subroutine computes the matrix difference 

C = A- B 


USAGE 

CALL SUBT(A,NA,B,NB,C,NC) 
Input Arguments 


Matrices: A,B 

Dimension arrays: NA,NB 

Output Arguments 

Matrices: C 

Dimension array: NC 


REMARK 

Matrices A and C may share the same storage space or matrices B and C may share the same 
storage space. 


7. MULT 
DESCRIPTION 

This subroutine computes the matrix product 

C = A B 
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USAGE 


CALL MULT(A,NA,B,NB,C,NC) 
Input Arguments 


Matrices: A,B 

Dimension arrays: NA,NB 

Output Arguments 

Matrix: C 

Dimension array: NC 


8. SCALE 
DESCRIPTION 

This subroutine multiplies every element of matrix A by S and stores the resulting value in B, 
that is, 



where S is a scalar. 

USAGE 

CALL SCALE(A,NA,B,NB,S) 

Input Arguments 

Matrix: A 

Dimension array: NA 

Scalar: S 

Note: If S is a constant, it must be written as a double precision constant (i.e., 2.D0, 0.D0, 

etc.). 

Output Arguments 

Matrix: B 

Dimension array: NB 

Note: A and B can be the same matrix. 
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9 . 


TRANP 


DESCRIPTION 

This subroutine rearranges the elements of matrix A so that 

B = A' 
or 

Bij = A ji 

USAGE 

CALL TRANP(A,NA,B,NB) 

Input Arguments 


Matrix: A 

Dimension array: NA 

Output Arguments 

Matrix : B 

Dimension array: NB 


10. INV 
DESCRIPTION 

This subroutine computes the matrix inverse of A and stores this inverse in A, that is, 


A = A ' 1 


Note that after the inversion is performed, the values stored in the original matrix A are destroyed 
and replaced by the corresponding elements of its inverse. 

USAGE 

CALL INV(A,NA,DET,DUM) 

Input Arguments 


Matrix: 

Dimension array: 


A 

NA 



Output Arguments 

Matrix: A, the inverse of the original A 

Scalar: DET, the determinant of A 

Dummy Argument 

Matrix : DUM, work vector of length 2*NA( 1 ) 

This subroutine is a slightly modified copy of the inverse routine given in the IBM scientific 
subroutine package. 

11. NORM 
DESCRIPTION 

This subroutine computes the norm of the matrix A as follows: 

|| A|| = min (max 2 Aj.-, max 
' j i J i 

USAGE 

CALL NORM(A,NA,ANORM) 

Input Arguments 

Matrix: A 

Dimension array: NA 

Output Arguments 

Scalar: ANORM 

12. UNITY 
DESCRIPTION 

This subroutine computes the unit matrix 

A = I 

USAGE 



CALL UNITY(A,NA) 



Input Argument 

Dimension array: NA 

Output Argument 

Matrix: A 

13. TRCE 
DESCRIPTION 

This subroutine computes the trace of the matrix A 

n 

TR= L a u 

i=i 

USAGE 

CALL TRCE(A,NA,TR) 

Input Arguments 

Matrix: A 

Dimension array: NA 

Output Argument 

Scalar: TR 

14. EQUATE 
DESCRIPTION 

This subroutine copies the values stored in matrix A into matrix B as follows: 

B = A 

USAGE 

CALL EQUATE(A,NA,B,NB) 

Input Arguments 

Matrix: A 

Dimension array: NA 
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Output Arguments 

Matrix: B 

Dimension array: NB 


15. JUXTC 
DESCRIPTION 

This subroutine takes the m X n matrix A, the m X p matrix B, and forms the m X (n+p) 
matrix 


C = [A B] 


USAGE 

CALL JUXTC(A,NA,B,NB,C,NC) 
Input Arguments 


Matrices: A,B 

Dimension arrays: NA,NB 

Output Arguments 

Matrix: C 

Dimension array: NC 


16. JUXTR 
DESCRIPTION 

This subroutine takes the m X n matrix A, the p X n matrix B, and forms the (m+p) X n 
matrix 


A 


C = 


B 


USAGE 

CALL JUXTR(A,NA,B,NB,C,NC) 
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Input Arguments 


Matrices: A,B 

Dimension arrays: NA,NB 

Output Arguments 

Matrix: C 

Dimension array: NC 


17. EAT 

DESCRIPTION 

This subroutine computes 


B = e At 


(*) 


and 


C = J e Ar dr 

o 

For a linear time-invariant system, the system equation is 


Then, 


or 


x = Ax + Gu 


x(t) = e At x 0 



drlGu 


x(t) = Bx q + CGu 


(*) 


See ASP manual, page 92, for reference. 

USAGE 

CALL EAT(A,NA,TT,B,NB,C,NC, DUMMY, KDUM) 

Input Arguments 

Matrix: A 

Dimension array: NA 

Scalar: TT 

where TT is the value of t used in equations * 
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I. 


Output Arguments 


Matrices: B,C 

Dimension arrays: NB,NC 

Dummy Arguments 

Matrix: DUMMY 

Constant: KDUM 


Note: KDUM contains the size of the DUMMY matrix, which must be at least 2*NA(1)*NA(2). 


18. ETPHI 
DESCRIPTION 

This subroutine computes the matrix exponential 

B = e At 


See ASP manual, page 92, and also EAT, page 18 of this manual for reference. 
USAGE 

CALL ETPHI(A,NA,TT,B,NB, DUMMY, KDUM) 

Input Arguments 


Matrix: 

A 

Dimension array: 

NA 

Scalar: 

TT 


where TT is the final value of time. 

Output Arguments 

Matrix : B 

Dimension array: NB 

Dummy Arguments 

Matrix: DUMMY 

Constant: KDUM 

Note: KDUM contains the size of the DUMMY matrix, which must be at least 2*NA(1)*NA(2). 
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19. AUG 


DESCRIPTION 

This subroutine computes 

and 


C = R’ 1 G' 


Z = 


-F 


+GR' 1 G' 


+H'QH +F' 

The matrices C and e are then used in RICAT to calculate the covariance and weighting matrices. 
These matrices arise from a linear system of the form 

x = Fx + Gu 


with output equation 


y = Hx 


and cost function 

3 = J (x'H'QHx + u'Ru)dt 
See ASP manual, page 2 1 2, for reference. 

In the special case where 


y = x 


then, 


-F 


Z = 


Q 


+GR* 1 G' 
+F' 


and the cost function is 

J = J* (x'Qx + u'Ru)dt 

A control index II is used to distinguish the two cases. 
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REMARKS 


The inputs to this program are the matrices F, G, RI, H, Q. 

(a) F must be square. 

(b) Q, R must be symmetric . 

(c) R must be invertible. 

The Fortran symbol for RT 1 is RI. 

USAGE 


CALL AUG(F,NF,G,NG,RI,NRI,H,NH,Q,NQ,C,NC,Z,NZ,II) 


Input Arguments 
Matrices: 

Dimension arrays: 
Control constant: 


Output Arguments 
Matrices: 

Dimension arrays: 


F,G,RI,H,Q 

NF,NG,NRI,NH,NQ 

II 

II ¥= 1 General case 

II = 1 Special case, H is not used in AUG 


C,Z 

NC,NZ 


20. RICAT 
DESCRIPTION 

This subroutine computes P(t) and K(t) by the following equations: 

P(t + r) = [6 2 1 +0 22 P(t)][0n +0 12 P(t)r 1 
K(t) = CP(t) 

See ASP manual, page 9, for reference. 

MOTIVATION 

A standard control problem will be used to illustrate how this matrix Riccati equation arises. 
Given the system equation. 


x = Fx + Gu 


the output equation, 
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y = Hx 


and the performance index, 

T 

J = J (x'H'QHx + u'Ru)dt + x'(T)H'S(T)Hx(T) 
o 

where Q,R,S are symmetric matrices and R is invertible. We wish to find a control law which 
minimizes the performance index J. Introducing the auxiliary variable X(t) into the system of 
equation, we have the following Euler-Lagrange equations, 4 

p -GR" 1 G' - 

\ -H'QH -F' 


pi 

= -z 

X 


_X_ 


which have for a solution 


x(t)“ 

= e Zt 

1 

o 

X 

1 

II 

I 

o 

X 

1 

_X(t)_ 


1 

0 

1 


— 1 

o 

1 


The optimal control law is 


u(t) = R' 1 G'X(t) 

Letting P(t) be a linear transformation from the state variable x(t) to the auxiliary variable X(t), 
that is. 


X(t) = P(t)x(t) 

we obtain from the Euler-Lagrange equation the following Riccati equation, 

-P = F'P + PF - PGR’ 1 G'P + H'QH 
where the initial condition for this differential equation is 

P(t) = H'S(T)H 


The optimal control, in terms of the state variable x(t), is 

u(t) = — K(t) x(t) = -R' 1 G'P(t)x(t) 

4 AUG computes Z rather than — Z, so that the exponentiation for d uses positive time increments. 
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and the optimal feedback gain K(t) is 


K(t) = — R" 1 G'P(t) 

Letting 

C = R M G' 

then, 

K(t) = CP(t) 

REMARKS 

1 . This subroutine will be terminated when 



or NC0NT(2) steps have been taken. 

2. Matrices P(t) and K(t) will be printed out every NCONT(l) steps, as controlled by 
NCONT(3). 

3. Matrices j, 0! 2 , 02 i > ^22 are submatrices of 0. Their dimensions are n X n where n is 
the order of the system (i.e., the dimension of the F matrix). They are partitioned from the 0 
matrix as follows: 



The Fortran symbol for 0 is PHI. 

USAGE 

CALL RICAT(PHI,NPHI,C,NC,NCONT,K,NK,PT,NPT,DUM,KDUM) 
Input Arguments 

Matrices: PHI, C, PT 

Dimension arrays: NPHI,NC,NPT 

Control array: NCONT(l) Number of steps per print 

NCONT(2) The maximum number of steps 

NCONT(3) Printout control 

1 -> no P, no K 

2 -»■ P only 

3 -*■ K only 

4 -»> P and K 
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Output Arguments 


Matrices: K,PT 

Dimension array: NK 

Dummy Arguments 

Matrix: DUM 

Constant: KDUM 


Note: KDUM contains the size of the DUMMY matrix which must be at least NPHI(1)**2. 

Note: PT is used for both input and output arguments. The initial value of P must be placed in 
PT before calling the subroutine. The value of P is updated every iteration in the subroutine until 
the final P is reached. This final P is one of the outputs of the subroutine. 


21. SAMPL 
DESCRIPTION 

Subroutine SAMPL calculates the covariance and weighting matrices associated with the 
discrete case of either the control problem or the filter problem. 

Consider the following filter problem. 

Given the system Xj +1 = 0Xj + u where u = gaussian random sequence 
with variance = Q, and observations yj = Hxj + v where v = gaussian 
random variable with variance = R. 


The optimum estimate of the state is (see p. 234 in the ASP manual) 

x i+1 ^xj + KiCyj-Hxj) 


where 


K i = 0P i H T (HP i H T + R) # 

P; + = Pj - PjH T (HPjH T + R) # HPj 

P i+ , = 0Pj + 0 T + Q 


# = pseudo inverse 


Here Pj is the solution of the matrix Ricatti equation, which is obtained by SAMPL. The subrou- 
tine has for inputs 0, H, Q, R, Pj, and for output, Pj +n and Kj +n _j where Pj +n is written over Pj. 

REMARKS 

1 . The routine will take n steps at a single call where n is an input parameter. Further, if P 
becomes constant, then the routine will stop and exit before completing the n steps. The actual 
test is as follows: 
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If 

2 . 

occurs. 


+ l)j|< 10’ 5 then 

The routine will print the value of Pj and/or Kj_, every j 
NCONT(3) controls which arrays are printed. 



exit. 

steps, and also when either exit 


USAGE 


CALL SAMPL(PHI,NPHI,H,NH,Q,NQ,R,NR,P,NP,K,NK,NCONT,DUM,KDUM) 


Input Arguments 


Matrices: 

Dimension arrays: 
Control arrays: 


Output Arguments 
Matrices: 

Dimension arrays: 

Dummy Arguments 

Matrix : 

Constant: 


PHI,H,Q,R,P 

NPHI ,NH,NQ,NR,NP 

NCONT 

NCONT( 1 ) = j = number of steps per print 
NCONT(2) = n = maximum number of steps 
NCONT(3) = print control 

1 no print 

2 print P only 

3 print K only 

4 print both P and K 


P,K 

NP,NK 


DUM 

KDUM 


Note: KDUM contains the size of the DUM matrix, which must be at least 6*NPHI(1)*NPHI(2). 


22. TRNSI 

This subroutine computes 

x(t) = e Ft x(t 0 ) + 

where 

u(t) = JR — Kx(t 0 + it 2 ) 
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and u is held constant for any interval specified by 

it 2 <t — t Q <(i+ l)t 2 i = 0,1,2,... 

The system output y(t) is given by 

y(t) = Hx(t) 

The state vector x and system outputs y are printed every tj intervals. Also t 2 must be a 
positive integral multiple of t t . The program terminates at t> tf. 

See ASP manual, pages 120-121, for reference. 

USAGE 

CALL TRNSI(F,NF,G,NG,J,NJ,R,NR,K,NK,H,NH,X,NX,T,DUMMY,KDUM) 

Input Arguments 

Matrices: F,G,J,R,K,H,X,T 

Dimension arrays: NF,NG,NJ,NR,NK,NH,NX 

Note: Dimension of T is 4 where 


Dummy Argument 

Matrix: DUMMY 

Constant: KDUM 

Note: KDUM contains the size of the dummy matrix, which must be at least 
4*NF(1)*NF(2). 

23. PSEUDO 
DESCRIPTION 

This subroutine computes the Moore-Penrose generalized inverse of the input matrix. It sets 
up a standard set of options for use by PSEU, which does the actual inversion. For details of the 
method, see PSEU, p.70. 


T 1 tt 

T2 t 2 

T3 t f 

T4 U 
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USAGE 


CALL PSEUDO(A,NA,B,NB,DUM,KDUM) 


Input Arguments 

Matrix: A 

Dimension array: NA 

Output Arguments 

Matrix: B = A# 

Dimension array: NB 

Dummy Arguments 

Matrix: DUM 

Constant: KDUM 


Note: KDUM contains the size of the dummy matrix, which must be at least 3*NA(1)*NA(2). 


24. DECGEN 
24a. DECSYM 

DESCRIPTION 

This subroutine decomposes a real matrix R with dimensions mXn and rank r<min(m,n) 
into two matrices H and G such that R = HG. Further, both H and G are of maximal rank, with 
dimensions m X r and r X n, respectively. It uses subroutine DECOM to provide matrices from 
which H and G can be computed. The writeup of DECOM, p. 85, describes the method in detail. 
Subroutine DECOM requires for input a matrix A which is positive semidefinite symmetric. Sub- 
routine DECGEN computes this matrix by letting A = RR^ or R^R, whichever is smaller, and uses 
the former if R is square. If the user knows that R is already positive semidefinite symmetric, 
this step may be omitted by a call to DECSYM, in which case A = R. 

USAGE 

CALL DECGEN(R,NR,G,NG,H,NH,DUM,KDUM) 
if R is general, or 

CALL DECSYM(R,NR,G,NG,H,NH,DUM,KDUM) 
if R is positive semidefinite symmetric. 
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Input Arguments 


Matrix: R 

Dimension array: NR 

Output Arguments 

Matrices: H,G 

Dimension arrays: NH,NG 

Dummy Arguments 

Matrix : DUM 

Constant: KDUM 


Note: KDUM contains the size of the DUM array, which must be at least 
7*min(NR( 1 ) 2 ,NR(2) 2 ). 


EXAMPLE USES OF VASP PROGRAM 


The examples given demonstrate directly the use of the principal subroutines EAT, ETPHI, 
AUG, RICAT, SAMPL, DECGEN, and PSEUDO. In addition, they exercise all of the subroutines 
except TRCE . They can be used to indicate whether the programs are working properly. They do 
not, however, provide an exhaustive test of the VASP program. 

The first example discusses the user’s main program in great detail to explain some of the 
system features. The remainder of the examples simply state the problem, and present the main 
program listing, the data listings, and the results. 


Example 1 - Transient Response 
A set of equations for a linear plant can be written as: 

x(t) = Fx (t) + Gu(t), x(0) = x Q 
y(t) = Hx(t) 

where x, u, and y are, respectively, the state, control, and observation vectors. The system, distribu- 
tion and observation matrices are F, G, and H, respectively. It is known that 

x(t) = e Ft x 0 + f e F ^ t ' r ^G(T)u(r)dT 

• o 

is the solution for x(t). If G and u are constant, then 
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x(t) = e Ft x 0 + J* e^*' 7 ") dr Gu 
o 

By letting s = t — r the integral becomes 

J e Fs d(— s) = J e Fs ds 

t o 

Thus, the solution to the system equation can be written 

x(t) = Bx q + CGu 
y(t) = Hx(t) 


where 


B = e Ft 


C = J e Fs ds 

o 

It is desired to generate the transient response of such a system in response to a given initial 
condition x Q and fix control u. In particular, given 



”l 

0 

0 


1 

0 

0 


F = 

0 

2 

0 

G = 

0 

1 

0 



_0 

0 

3_ 


_0 

0 

1_ 
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1 



1 


1 

1 

1 






H = 



u - 

0 

x o 

= 

2 


0 

1 

0 


_0_ 



_3_ 


find x(t) for 0< t < 2.0. Also print x(t) and y(t) every 0.01 second. 

The user’s main program to solve this problem is shown in figure 1(a), the corresponding data 
deck is shown in figure 1(b), where each line represents one card, and the beginning of the results is 
presented in figure 1(c). 
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OJ 

o 


0001 

DIMENSION F(3,3) ,NF(2) ,G(3,3),NG(2) ,H(2,3> ,NH(2 ) ,B(3,3) , 
XNB(2) ,C(3,3),NC(2),XO(3) , NX0 ( 2 ) , XT ( 3 ) , NXT ( 2 ) , VI ( 3 ) , NV 1 ( 2 ) 

yu 9 m.uuon\.Ain. 7 » w»im iiin.Mii/'n vt/o\ mvt n i unot 

f 

000 2 

0003 

0004 

DOUBLE PRECISION F , G ,H ,B , C ,X0, U, XT , YT , TT , DELTAT ,TFINAL , VI 

COMMON /MAX/ MAXRC 

COMMON /LINES/NLP • L IN-.TT TL F( 73 ) 

,V2,A1,W 

0005 

0006 

0007 

MAXRC=9 
CALL RDTITL 


000 8 

0009 

0010 

READ (5,100) TT,0ELTAT ,TF INAL 
100 FORMAT (8F10.2) 

CALL READ ( 5 • F • NF . G • NG . H » NH . IJ • NU. XO .NXO ) 


0011 

101 FORMAT! 1H0 59X, ‘TIME RESPONSE', /• T IME • ,22X ,• STATE • ,54X, 
1 'OUTPUT*, /• TT'»14X,'XT(1)',11X,'XT(2)'»11X,'XT(3)'»21X, 

1 11 V.IVTnil.llY JVTI 1 LIL 

' YT ( 1 ) • , 

0012 

0013 

0014 

102 FORMAT ( 1H0 F10 . 2 ,6X , 3E 16 . 7 , 10X, 3E1 6.7 ) 
CALL LNCNT (4) 

WRITE ( 6 . 101 ) 


0015 

0016 
on l 7 

n x- t -t 1 v 7 iv t f 

10 CALL EAT ( F , NF, TT,B,NB,C,NC , W, 18 ) 

CALL MULT( B , NB , XO ,NX0, VI ,NV1 ) 

r a_j__ j_ kii ii Tir kir* r Kit' a i ma i i _ _ 


0018 

0019 

0020 

wMLl. H UL I ILf NU y \j y Nb 7 M 1 y Iih i I 

CALL MULT(A1,MA1,U,NU,V2,NV2) 
CALL AD0( VI , NV1 , V2 ,NV2 , XT, NXT ) 

CAM Mill TlH.MH.yT.MVT.VT.NVT1 


0021 
0022 
on 0 3 

un u t_ t t n 7 1 \'D 7 A t 7 ItA t 7 T I 7*rrr I f 

CALL LNCNT (2) 

WRITE (6, 102) TT, ( XT ( I) ,1=1,3), ( YTC I ) ,1=1,2) 


0024 

0025 
on?A 

IF(TT.GT.TFINAL) STOP 
GO TO 10 

PMn 


(a) User’s main program. 


Figure 1.— Example 1. 



TEST 

PROGRAM 1 


GENERATES TRANSIENT USING EAT 

.01 

.01 


2.00 

F 

3 

3 


1.0 

0.0 


0.0 

0.0 

2.0 


0.0 

0.0 

0.0 


3.0 

G 

3 

3 


1.0 

0.0 


0.0 

0.0 

1.0 


0.0 

0.0 

0.0 


1.0 

H 

2 

3 


1.0 

1.0 


1.0 

0.0 

1.0 


0.0 

u 

3 

1 


1.0 




0.0 




0.0 




xo 

3 

1 


1.0 




2 .0 




3.0 





(b) Data deck. 


Figure 1 Continued. 



& TEST PROGRAM 1 GENERATES TRANSIENT USING EAT 


VASP PROGRAM 


F MATRIX 3 ROWS 3 COLUMNS 

1 , 00000000 00 0^-0 9 - r 0 

0.0 2.00000000 00 0.0 

0.0 0.0 3 • 0000000 D 00 

G MATRIX 3 ROWS 3 COLUMNS 

1 • OOOOOOOD 00 0.0 0.0 

■O^Q 1 . OOOOOOOD 00 0^9 

0.0 0.0 1.00000000 00 

H MATRIX ... 2 ROfeLS _ 3 COLUMNS 

1. OOOOOOOD 00 1. OOOOOOOD 00 1. OOOOOOOD 00 

0.0 1. OOOOOOOD 00 0.0 


U MATRIX 3 ROWS 1 COLUMNS 

1. OOOOOOOD 00 

....CL. a 

0.0 

XD MATR-LX 3-ROWS 1 COLUMNS 

1. OOOOOOOD 00 
2. OOOOOOOD 00 
3 . OOOOOOOD 00 


T |MP 


5IATC. 



TIME RESPONSE 

niiTmiT 


TT 


XT ( 1 ) 

XT ( 2 ) 

XT ( 3 ) 


YT ( 1 ) 

YT ( 2 > YT ( 3 ) 

— 

cuai . 

_ ... 0*10201000 

01 

0*20404030. 01. 

0*309 1 26 4D 

i 

01 

o*6i5ia67.a 

-01 - 

0,20404030-04 


0.02 

0 .10404030 

01 

0 • 208 1622D 

01 

0.31R5510D 

01 

0 • 63075340 

01 

0 • 208 1622D 01 


o 

• 

o 

0 • 1060909D 

01 

0.21 23673D 

01 

0 . 32 82 52 3D 

01 

0 . 64671050 

01 

0.21236730 01 


_a.Q4 

0,1 Oft 1 6 77 n 

01 

0,21665740 

01 

0.3382491D 

U1 

Q ♦ 6.6 20-6-1160 

OL ... 

0*2 166-53 40- -Cll 


0.05 

0.11025420 

01 

0.2210342 D 

01 

0 .34855030 

01 

0. 679 83870 

01 

0.22 10342D 01 


0.06 

0.1123673D 

01 

0 • 22 54994D 

01 

0.3591652D 

01 

0 • 69703 1 9D 

01 

0. 22549940 01 


0,07 

0.H450L60 

01 

0.23005480 

01 

0.370 10340 

01 

0 • 7 146598D 

01 

0.23.005480 0.1 _ 







(c) Output 






Figure 1.— Concluded. 
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The user’s main program— This program will be discussed, statement by statement, using the 

line numbers on figure 1(a) as a reference. 

Lines 1 and 2. These two statements allocate the necessary storage for the variables to be used and 
define them as double precision. Also, the dimension arrays NF, NG, etc., are allocated storage. 
The dimensionality of F, G, etc., could have been included in the double precision statement 
instead of the dimension statement, and they could have been dimensioned as F(9) instead of 
F(3,3). The W array has been set up for dummy storage, and is dimensioned 18, as required 
by the EAT subroutine. 

Lines 3 and 4. Common variables to be needed later are made available to the program. Although 
the variables listed in line 4 are not needed in this program, they are shown for reference. 

Line 5. Since the basic matrices are (3,3), MAXRC is set to 9, to prevent overfilling the matrices. 
Note this will not protect from overfilling the arrays XO, XT, etc., since they are expected to 
be 3 X 1 vectors, and are dimensioned 3. 

Line 6. This statement reads the first card of the data deck (see fig. 1(b)), places its contents in the 
TITLE array, and prints the first line of the output (see fig. 1(c)). 

Lines 8 and 9. The initial time, the time increment, and final time are read from the second data 
card. 

Line 10. The arrays F, G, H, U, and XO are read from the remainder of the data deck, and are 

printed (fig. 1(c)). Note that the dimensions used by the program are those given on the header 
card for each matrix. If these were specified as (2,2) this same main program would solve a 
second-order problem, rather than the third-order problem. 

If the initial conditions were already stored in the XO array and you did not wish to disturb 
them, then the header card for the XO array would contain only the matrix title, no dimen- 
sions, and the associated data cards would be omitted. The matrix XO would still be printed. 

Line 1 1. Line 1 1 contains the information to head the main output. 

Line 1 2. Line 12 is the data format. For this program the transient output was printed using the 
programmers write statement rather than PRNT. The use of PRNT for this purpose is shown 
in the third example, p. 40. 

Line 13. Line 13 tells the line counter that the program will print 4 lines. 

Line 14. Line 14 does the actual printing. 

Lines 15 through 25. Lines 15 through 25 form a loop which increments TT (line 23) and stops 
when TT is large enough (line 24). 

Line 15. Line 15 computes the B and C matrices for time TT. When C is computed, the limits of 
the integral are 0 and the present TT. Note that W is specified for dummy storage and the 
“18” tells EAT the size of W. 
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Line 16. Line 1 6 computes BXO and stores the result in VI. Array VI is set up for the 

programmers working storage. Since W is also available at this point in the program, it could 
have been used instead of V 1 if desired. 

Line 17. Line 17 computes CG and stores the result in Al, another working storage array. 

Line 18. Line 18 computes (CG)U and stores the result in V2, still another working storage array. 
Note that MULT obtains the product CG from A 1. 

Line 19. Line 19 adds VI and V2 to obtain XT. Since the ADD subroutine allows the matrices to 
be repeated in the call, the array VI could have been eliminated, then line 16 would have 
stored its results in XT. Line 19, then, would have added XT and V2 to obtain the complete 
XT. 

Line 20. Line 20 multiplies H times X to obtain Y. 

Line 21 . Line 21 tells the counter we are going to print 2 lines. If this will not fit on the present 

page, LNCNT will advance to the next page, print the title as on the first line of the first page of 
output, and increment the line counter to allow for the paging and the two lines about to be 
printed. 

Line 22. Line 22 prints the variables XT and YT, skipping a line between each print line, as 
required by the 1HO in FORMAT 102. Note that YT(3) is not printed. 


Example 2 — Transient Response Using TRNSI 

This example uses the same equations as Example 1 , except that u is piecewise constant, 
that is, 


u(t) = JR - Kx(t 0 + it 2 ) it 2 < t - t Q < (i + l)t 2 

where i is a non-negative integer and J, R, K are constant matrices. The first term, JR, represents 
a forcing function and the second, Kx, is a feedback term. (See ASP manual, p. 121, for detailed 
explanation.) 

It is desired to generate the transient response of such a system in response to a given initial 
condition x 0 and a time varying control u. In particular, given 



"0 

0 

0 

0 

0 “ 


1 

0 - 


0 

5 

0 

0 

0 


1 

1 

F = 

0 

0 

1 

1 

0 

G = 

0 

0 


0 

0 

0 

1 

0 


1 

0 


_0 

0 

0 

0 

2_ 


0 

1 

J = 

(-01 

LoJ 


R = 

[0] 

K 

-p 0 

|_0 3 


0.5 

0 
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H = 


1 

0 

0 

0 

0 

1 

0 


0 0 
1 0 
0 1 
0 0 
0 0 
0 0 
1 0 


0 0 
0 0 
0 0 
1 0 
0 1 
0 0 
1 0 


1 

1 

1 

0 

1 


tj = 0.5 sec 
t 2 = 2 se£ 
tf = 3.5 tec 
t Q = 0 sec 


The system is monitored at intervals tj , while the control u(t) is changed only at sampling 
intervals t 2 (t 2 must be a positive integral multiple of tj ). Specifically, the control u(t) is updated 
by the equation: 


u(t) = JR — Kx(t 0 + it 2 ) it 2 < t — t Q < (i + 1 )t 2 

The x,y vectors are computed at time intervals tj , and these vectors together with the time t, and 
the control u (for the subsequent time interval) are printed out each time. The problem terminates 
when the final time tf is reached. The matrix T has elements t 2 , t 2 , tf, t Q in that order. 

The user’s main program to solve this problem is shown in figure 2(a), the corresponding data 
deck is shown in figure 2(b), and the results are presented in figure 2(c). 


100 DIMENSION F(5t5)»NF(2)*G(5*2)*NG(2)*J(?»l) »N J ( 2 ) » R ( I * 1 ) »WRJ2 ]_,_K ( ?t 

200 15 ) * NK ( 2 > * H ( 7 » 5 ) » NH ( 2 ) * X ( 5 ) *N X ( 2 ) * T ( 4 > ,NT<2) , DUMMY! 100) 

300 DOUBLE PRECISION F ,G , J , R t K , H ,X , T , DUMMY 

i 00 COMMON/ FORM/MEPR t FM T1 ( 6) ,FMT2( 6) 

500 CALL RDTITL 

600 CALL READ(5,F,NF,GtNG,J,NJtR»NR,K,NK) 

700 CALL READ(3,H,NH,X,NXfT,NT) 

BOO CALL TRNSi'(F,NF f G ? MGtJtNj,R,MR,K,NK f H t NH,x'VNX/f,ni7MMY,100) 

900 STOP 

1000 .END ___ _ .. . 


(a) User’s main program. 


Figure 2.— Example 2. 
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TEST 

PROGRAM 


GENERATES 

TRANSIENT 

USING 

TRNSI 

F 

5 

5 





0, 

0. 


0. 

0. 

0. 


0. 

.5 


0. 

0. 

0. 


0. 

0. 


1. 

1. 

0. 


0. 

0. 


0. 

1 . 

0. 


0, 

0. 


0. 

0. 

2. 


G 

5 

2 





1. 

0. 






i. 

1. 






0. 

0. 






l . 

0. 






0. 

1 . 






j 

2 

1 





0. 







0. 







R 

1 

1 





0. 







K 

2 

5 





1 . 

0. 


.5 

0. 

2. 


0. 

3. 


0. 

1. 

0. 


H 

7 

~ 5 ~ 





1. 

0. 


0. 

0. 

0. 


0. 

1. 


0. 

0. 

0. 


0. 

0. 


1. 

0. 

0. 


0. 

0. 


0. 

1. 

0. 


0. 

0. 


0. 

0. 

1. 


1. 

0. 


0. 

0. 

0. 


0. 

1. 


0. 

1. 

0. 


X 

5 

1 





1 . 







1 . 







1 . 







0. 







1 . 







T 

4 

1 





.5 







2 . 







3.5 








0 . 


(b) Data deck. 
Figure 2.— Continued. 


36 



TEST PROGRAM 


GENERATES TRANSIENT USING TRNSI 


VASP PROGRAM 


F 

MATRIX 5 ROWS 

5 COLUMNS 




0.0 

0.0 

0.0 

0.0 


670 ” 

0.0 

5 . OOOOOOOD-O 1 

0.0 

o 

• 

o 


o 

• 

o 

0.0 

0.0 

I • OOOOOOOD 00 

1. OOOOOOOD 

00 

0.0 

6.0 

0.0 

0.0 

1. OOOOOOOD 

oo 

0.0 

0.0 

0.0 

0.0 

o 

« 

o 


2. OOOOOOOD 00 


G MA T R I X 5 R O WS 2 COLUMN ' S 


1. OOOOOOOD 00 
1. OOOOOOOD 00 

0.0 

1. OOOOOOOD 00 







o.o 

1. OOOOOOOD 00 
0.0 

0.0 

0.0 

1. OOOOOOOD 00 







J 

0.0 

MATRIX 


2 ROWS 

1 

COLUMNS 





0.0 










R 

MATRIX 


1 ROWS 

1 

COLUMNS 





o 

• 

o 










K 

MATRIX 


2 ROWS 

5 

COLUMNS 





1. OOOOOOOD 00 
0.0 

0.0 

3. OOOOOOOD 00 

5. 0000000D-01 
0.0 

0.0 

1. OOOOOOOD 

00 

2. OOOOOOOD 
0.0 

00 

H MATRIX 

1. OOOOOOOD 00 
0.0 

7 ROWS 

0.0 

1. OOOOOOOD 00 

r 

0.0 

0.0 

COLUMNS 

0.0 

0.0 


0.0 

0.0 


0.0 

0.0 

0.0 


0.0 

0.0 

0.0 


1. OOOOOOOD 00 
0.0 
0.0 

0.0 

1. OOOOOOOD 
0.0 

00 

0.0 

0.0 

1. OOOOOOOD 

00 

1. OOOOOOOD 00 
0.0 

0.0 

I. OOOOOOOD 00 

0.0 
0 . 0 


0.0 

1. OOOOOOOD 

00 

0.0 

0.0 



u> 

-j 


(c) Output 

Figure 2.- Continued. 



u> 

00 


VASP PROGRAM 


TEST PROGRAM GENERATES TRANSIENT USING TRNSI 

X MA T R IX 5 -RO Vf S r-C O t U MN S 

1 . OOOOOOOD 00 

l.OOOOOOOP 00 

1. OOOOOOOD 00 

0.0 

1. OOOOOOOD 00 


T MATRIX 4 ROWS 1 COLUMNS 

5 . OOOOOOOD- 01 
2. OOOOOOOD 00 
3 • 5000000D 00 
0.0 


E 

0.0 

0.0 

MATRIX 

5 ROWS 

0.0 

5 . OOOOOOOD-O 1 

5 ' COLUMNS 

0.0 

0.0 

o o 

• • 

o o 

o o 

• • 

o o 

OTO 


TT70 

1 .TKTOOO'OOl) “00 ' 

1. 00000000 00 

0.0 

0.0 


O 

• 

o 

O 

• 

c 

1. OOOOOOOD 00 

0.0 

0.0 


0.0 

o 

• 

o 

0.0 

2. OOOOOOOD 00 

EAT 

MATRIX 

5 ROWS 

5 COLUMNS 



1. OOOOOOOD 00 

0.0 

0.0 

0.0 

0.0 

0.0 


1 • 2840254D 00 

0.0 

0.0 

0.0 

0.0 


0.0 

1.648721 30 00 

8.2436069D-01 

0.0 

0.0 


0.0 

0.0 

1.648721 3D 00 

0.0 

O 

• 

o 


0.0 

0.0 

0.0 

2.7182819D 00 

INT 

MATRIX 

5 ROWS 

5 COLUMNS 



5 . 0000002 D-01 

0.0 

o 

• 

o 

0.0 

O 

• 

O 

0.0 


5 .68050860-01 

0.0 

0.0 

0.0 

0.0 


0.0 

6.48721310-01 

1.7563938D-01 

0.0 

0.0 


0.0 

0.0 

6. 48721 3 1 D-0 1 

0.0 

O 

• 

o 


0.0 

0.0 

0.0 

8.59140970-01 


(c) Output - Continued. 
Figure 2.— Continued. 



TEST PROGRAM GENERATES TRANSIENT USING TRNSI 




VASP PROGRAM 




TRANSIENT 

RESPONSE, * 

INDICATES CONTROL CHANGES 










TIME 

FIRST 5 

ELEMENTS CONTAIN X, NEXT 7 ELEMENTS CONTAIN Y 

= HX, LAST 

2 

ELEMENTS CONTAIN U =JR -KX 


o 

• 

o 

* 

1 • OOOOOOOD 

00 

1.00000000 

00 1.00000000 

00 

0.0 


1. OOOOOOOD 

00 

1. OOOOOOOD 

00 

1.00000000 

00 


X . OOOOOOOD 

00 

“0715 

1. OOOOOOOD 

ucr 

1 .OOOOOOOD 

uir 

1. OOOOOOOD 

W 

-3. 500000075" 

ui r 

-3. OOOOOOOD 

W 

0.50 

-7.5000008D-01 

-2.40830520 

00 1.03398350 

00 

-2.2705246D 

00 

1.4085903D- 

-01 

-7 .50000080- 

-01 

-2.40830520 

00 


1 .03398350 00 

-2.27052460 

Uo 1 .408 5903D- 


-7 • 5060OO 8D- 


-4.67682970 

"W 

-3.5000000D 

THT 

-3. ooOOoooD 

"TO* 

I. 00 

-2 • 5000002 D 

00 

-6.78465570 

00 -7.8171 846D- 

-01 

-6 .0 139868D 

00 

-2.1945284D 

00 

-2.50000020 

00 

-6.78465570 

00 


-7.81718460- 

-01 

-6.01 39868D 

00 -2.1 9452840 

00 

-2 .50000020 

00 

-1 . 2798642D 

01 

-3 .50000000 

w 

-3. OOOOOOOD 

00 

1.50 

-4 .2 500002D 

00 

-1. 2404001D 

01 -6.8612680D 

00 

-1 .218591 3D 

01 

-8 . 5427698D 

00 

-4. 2500002D 

00 

-1.24040010 

01 


-6.861 2680D 

00 

-1.21859130 

01 -8.54276980 

00 

-4.2500002D 

oTT 

-2.4589914D 

TH7 

-3.5000000D 

“W 

-3. OOOOOOOD 

THT 

* 2.00 

-6.0000003D 

00 

-1.9619383D 

01 -2.197 2644D 

01 

-2.2361699D 

01 

-2.57990800 

01 

-6.00000030 

00 

-1.96193830 

01 


-2.19726440 

01 

-2.2361699D 

01 -2 . 5 799080D 

01 

-6.0000003D 

00 

-4.19810820 

01 

6.8 584482D 

01 

8.1219849D 

01 

2.50 

2.82922420 

01 

5 • 9904692D 

01 -4.2614736D 

01 

7.6240057D 

00 

-3.49072850-01 

2.82922420 

01 

5.9904692D 

01 


-4.26147360 

01 

7.62400570 

00 -3.49872850- 

-01 

2 .8292242 D 

01 

6.75286970 

01 

6.85844820 

01 

8.12198490 

01 

. 

o 

o 

6.2584484D 

01 

1.62015630 

02 -5.19 2 8756D 

01 

5 • 70620750 

01 

6.88282470 

01 

6.25844840 

01 

1.62015630 

02 


-5 • 1928756D 

01 

5.7 06207 5D 

01 6.882 R247D 

01 

6.25844840 

01 

2.19077700 

02 

6.85844820 

01 

8 • 12 19849D 

01 

3.50 

9.68767270 

01 

2.9312866D 

02 -2.6530179D 

01 

1.3857167D 

02 

2 • 5687388D 

02 

9. 6876727D 

01 

2.93128660 

02 


-2.65301790 

01 

1.38571670 

02 2.5 687 3 88D 

02 

9.68767270 

01 

4. 31 70034D 

02 

6.85844820 

01 

8 • 12 19849D 01 


co 

'•O 


(c) Output — Concluded, 
Figure 2.- Concluded. 



The user’s main program — A brief explanation of the statements using line numbers on 
figure 2(a) as reference follows: 

Lines 1, 2, and 3. Lines 1, 2, and 3 allocate storage, same as lines 1 and 2 of example 1. 

Line 4. Common variables to be needed later are made available to the program. 

Line 5. This statement reads the first card of the data deck (see fig. 2(b)), places its contents in 
TITLE array and prints the first line of the output (see fig. 2(c)). 

Lines 6 and 7. The matrices F, G, J, R, K, H, X, and T are read in from data deck (see fig. 2(b)) 
and are printed. 

Line 8. Line 8 calls the TRNSI subprogram, performs the computation, and prints outputs as 
explained in the example. 


Example 3 — An Optimum Control Problem 


Given a system 


x = Fx + Gu y = Hx x(o) = x Q 

where x, u, and y are, respectively, the state, control, and observation vectors. The system, 
distribution, and observation matrices are F, G, and H, respectively. 

We wish to define an optimal control u(t), where u(t) = -Kx(t), so as to minimize the 
performance index 

J= J* (x'H'QHx + U'RU)dt 


The solution to this problem is 


K = R' 1 G'P 


where P is the solution of the matrix Ricatti equation. 

The VASP program finds P by means of the subroutines AUG, ETPHI, and RICAT, as follows. 
First, subroutine AUG is used to generate the matrices 


Z = 


-F 

H'QH 


GR' 1 G' 
F' 


and 


C = R _1 G' 


(Note: This is the negative of the Z given on page 212 of the ASP manual.) Subroutine ETPHI is 
then used to compute the special transition matrix 
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e = 


Q\l 0\ 2 

0 2 1 $22 


= e ZT 


Finally, the P matrix is computed by subroutine RICAT for a time increment of r, by repeated 
application of the formula 


P(t + r) - [02 1 + 02 2^(0] [011 + 01 2^(01 1 

The computation is repeated for several steps, until P(t + r) « P(t), which is the desired 
solution. Subroutine RICAT will also stop after a specified number of steps, if P has not con- 
verged to a solution. Finally, having P and K, we can compute the transient response of the system 
with optimum feedback from any desired initial condition. The differential equation becomes 

x = Fx — GKx = (F — GK)x = F*x 


and the solution is 


x(t) = e^*^x 


The time history of the control is 


u(t) = — Kx(t) 


An alternate solution, used in this example, is to first calculate the transition matrix 


A2 = e 


F*t 


i 


where is the time increment at which the solution is desired, then compute 

x(t + r,) = A2 x(t), x(0) = x o 

The listing of a main program to solve this problem is given in figure 3(a), the data for a particular 
case is given in figure 3(b), and the first part of the results is given in figure 3(c). In this problem, 
H = I so the special case of AUG is used. As a result, H is not used in AUG, and need not have 
been used as an input. 
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to 


0001 


0002 


0003 

000 4 


DIMENSION NF ( 2 ) , NG ( 2 ) , NFS TAR ( 2 ) , NC ( 2) ,NCK(2 ) , NO { 2 ) , NR ( 2 ) , NR I ( 2 ) , 

1 NZ I 2 ) » N P H I ( 2) ,NP0(2) ,MPT(2) ,NXO( 2) , NXT ( 2 ) , N A 1 ( 2 ) , NA2 ( 2) ,NH(2) , 

2 N0E L TC( 2 ) , NXX(2) , Nr , nNT(3) ,lAaL 24 .,- L 4-64 


DOUBLE PRECISION 0 ( 3 , 3) , R ( 1 , 1 ) , R I ( 1 , 1 ) » F (3 ,3 ) ,G (3 , 1 ) , FSTAR ( 3 , 3 ) t 

1 C( 1,3) ,CK< 1,3) ,Z (6,6) ,PHI(6,6) , P0<3,3) ,PT(3,3) ,XO(3) ,XT(3) ,XX(3) , 

2 DET 

COMMON / MAX / MAXRC 


0005 

0006 

0007 

0008 

0009 

0010 


0011 

0012 

0013 

0014 

0 0 1 5 


0016 

0017 

0018 

0019 

0020 
0021 
0022 
0023 
0024- 

0025 

0026 
0027 


CALL RDTITL 
MA XRC = 3 6 

KPII M^ 7 ? — - - 

READ ( 5, 100) TT,DELT1,DELT2,TFINAL 

CALL READ ( 5 ,F ,NF , G, NG ,X0 ,MXO, H,NH,0,N0 ) 

GALL— READ < 2 , R , MR , PO , NP -Q , ft ,A1R-, R^- NR r R -, N R ) 


100 FORMA T ( 8F 10 « 2 ) 

101 FORMAT I 1H 59X , 'TIME RESPONSE • /5X, • T IME • , 22 X , • STAT E • ,43X , 

1 »OUTPUT» /5X i--XT»-r l 4 X , « X TWMlXr^ZimiXt 1 XT-4-3-4 1 -r24-Xr- t O&LTG 1 )- 

102 FORMAT ( 1H F 10 . 2 , 6X , 3E 16 . 7 , 10X , 2 El 6 . 7 , 10 X , F 16 . 7 ) 

NCONT ( 1 ) = 1 

MG & NT < 2 ) = 1 01 NT ( TF INAL/PEL T T-L+ 4 


NCONT ( 3 ) =4 

CALL EQUATE (XO,NXO»XT»NXT) 

CALI EQUATELPQ-, NPa^PI+MPX ) 

CALL EQUATE I R,NR ,RI ,NR I ) 
CALL I NV ( R I , NR I , DET, L ) 

CA1 I ALLG-t E,NF,G,NG+B-I-,AI 
CALL ETPHI (Z,NZ,DELT1 ,PHI 
CALL RICATIPHI " k ' 


0028 


ETPHI (Z,NZ,DELT1 ,PHI ,NPHI , DUMMY, KOUM) 

RICATIPHI , NPH I ,C,NC , NCONT, CK,NCK ,PT,NPT, DUMMY, KOI IM) 
UMLL MULTI G , NG, CK,NCK , A 1,NA1 ) 

CALL SCALE! A 1 , NA 1 , A2 ,NA2 , -1. DO) 

CALL ADD(F,NF,A2,NA2, FSTAR , NFS TAR ) 


r ai i ctdu 




(a) User’s main program. 
Figure 3.— Example 3. 



0029 

-miQ 

0031 
00 32 

-00*3 

0034 

0035 

-0034 

0037 

0038 

0039 

,0040 

0041 

0042 

004 ? 

0044 

0045 

004 4 

0047 

0048 


CALL LNCNT ( 2 ) 

— -W4UT-E- ( 4t1 0 6) - 

106 FORMAT CO TRANSITION MATRIX') 

CALL PRMT ( A2, NA2 , ' A2',l) 

XAL-U-U0 C -N -T4 1OO ) 

CALL LNCNT ( 3 ) 

WR ITE ( 6 » 101 ) 

-^ALlr-MyL TICK t NCK yXTtAFXTrOF^TC-rN&^LTOf 

DELTC! 1)=-1.0*DELTC< 1) 

CALL LNCNT ( 1 ) 

WR I TE ( 6 » 102 ) TT , XT, DELTC 

200 CAM FOIIATEUT.»NXT,XX,NXX) 

CALL MULTI A2,NA2, XX, NXX,XT,NXT) 

CALL MULT(CK,NCK,XT,NXT, DELTC , NDELTC ) 

TT = TT+DFLT2 

WRI TE I 6» 102 ) TT, XT, DELTC 

I TT . GE.T F I N AL) 5WP 

GO TO 200 
END 


u> 


(a) User’s main program - Concluded. 
Figure 3.— Continued. 



TEST PROGRAM 2A 


GENERATES TRANSIENT USING X ( 1+1 ) =E XP ( F*T) *X ( I ) 

0.0 

1.0 


0.01 3.5 

F 

3 

3 


-0.2767 

1.0 


-0.0372 

-17.0872 

-0 .1785 

-12.1983 

0.0 

0.0 


-6.67 

G 

3 

1 


0.0 




0.0 




6.67 




xo 

3 

1 


1.0 




0.0 




0.0 




H 

3 

3 


1 .0 

0.0 


0.0 

0 .0 

1.0 


0.0 

0.0 

0.0 


1.0 

0 

3 

3 


0.2 

0.0 


0.0 

0 .0 

0.2 


0.0 

0.0 

0.0 


6.0 

R 

1 

1 


1 .0 




PO 

3 

3 


0.0 

0.0 


0.0 

0.0 

0.0 


0.0 

0.0 

0.0 


O 

• 

o 


(b) Data deck. 


Figure 3.— Continued. 





TFST PROGRAM 2A. GENERATES TRANSIENT USING X ( I + 1 ) = F XP ( F *T ) *X ( I ) VASP PROGRAM 

F MATRIX 3 ROWS 3 COLUMNS 

_-L2„.76 700000^01 4 .OGOOGOO© 00 - -3.7 2000000-02 

-1 .70872000 01 -1 .78500000-01 - 1 . 2 198300D 01 

0.0 0.0 -6.67000000 00 

G MATRIX 3 ROWS 1 COLUMNS 

0.0 

6.67000000 00 

. XO MATRIX 3 ROWS 1 COLUMNS 

1 . OOOOOOOD 00 

0.0 

--cua 

H MATRIX 3 ROWS 3 COLUMNS 


1.00-000000 0© 

0.0 

0.0 


0.0 

1. OOOOOOOD 00 

0.0 


0.0 

0.0 

1. OOOOOOOD 00 

0 MATRIX 

3 ROWS 

3 

COLUMNS 

2 . OOOOOOOD-Ol 

0.0 

0.0 


0.0 

2. OOOOOOOD-Ol 

0.0 


0.0 

0.0 

0.0 


R . . . MAJR.IX 

1. OOOOOOOD 00 

1 ROWS 

1 

COLUMNS 

PO MATRIX 

3 ROWS 

3 

COLUMNS 

0.0 

0.0 

0.0 


0.0 

0.0 

0.0 


0 . 0 - 

0.0 

0.0 



1 ITERATIONS 


Lr, 


(c) Output. 

Figure 3.- Continued. 



GENERATES TRANSIENT USING X ( I + 1 ) = E XP ( F *T ) *X ( I ) 


VASP PROGRAM 


& TEST PROGRAM 2A 

K ( T ) MATRIX 1 ROWS 

7.36037320-01 -3 .63599430-0 1 

P.tJl MATRIX 3- -ROWS 

6.67994270-01 -2.16153470-02 

-2.16153470-02 5 . 5929 130D-02 

1.10350420-01 -5.45126590-02 

2 ITERATIONS 

K ( T ) MATRIX 1 ROWS 

— 7.S4633-140--O4 -^3*691 79860=0 1- 

PIT) MATRIX 3 ROWS 

6*76369130-01 -2.17721540-02 

-2.17721540-02 5.6466312D-02 

1.13138400-01 -5.53493070-02 

■ 3 -ITERATIONS 

KIT) MATRIX 1 ROWS 

7 *5462 2 8 30-01 -3*69 203940-01 

PIT) MATRIX 3 ROWS 

— 0* 7 6403960-01 2*176439 10 - 02 - 

-2.17641910-02 5.64703680-02 

1 • 1 3 1 3686D-01 -5.53529140-02 

4 ITERATIONS 

KIT) MATRIX 1 ROWS 

...;4r-54622 750-04 -3. 69604180-01 


3 COLUMNS 
5.16454380-01 

3 COLUMNS 
1.10350420-01 
5. 4512659 D- 02 
7.74294420-02 


3 COLUMNS 
5.30329180-01 

3 COLUMNS 
1.13138400-01 
5.53493070-02 
7 .95096220-02 


3 COLUMNS 
5. 303612GD— 01 

3 COLUMNS 
1*13136860-01 
-5.53529140-02 
7 .95144230-02 


3 COLUMNS 
5.30361530-01 


(c) Output - Continued. 
Figure 3.- Continued. 



P ( T ) MATRIX 3 ROWS 

6-.7 640 4 04D -Q 1 -2,1 7441 710-92 

-2.17641710-02 5 . 6470397D-02 

1 . 131 3684D-0 1 -5. 53529500-02 


| FSTR MATRIX 3 ROWS 

I -2.7670000D-01 1.00000000 00 

, ^L.7Q&7200Q Ql -U7B-500000-ai 

-5.0333338D 00 2.4625919D 00 

TRANSITION MATRI X 

A2 MATRIX 3 ROWS 

q QA 4 Q 41 2 0— 0 1 9*9651 89! D-03 

-1.67390330-01 9. 95924740-01 

-4 . 9775 05 5 D- 02 2.31282830-02 


(c) Output - Continued. 
Figure 3.- Continued. 


3 COLUMNS 
1,131 34840-01 
■5.53529500-02 
7.9514472D-02 

3 COLUMNS 
■3. 7200000 D-02 
- 1^219 8 300D 01 
■1.02075110 01 


3 COLUMNS 
9-^41 4194 80-04 
1.15736830-01 
9.0157826D-01 



VASP PROGRAM 


TEST PROGRAM 2A GENERATES TRANSIENT USING X ( 1+ 1 ) =EXP ( F*T ) *X ( I ) 


00 


TIME RESPONSE 


IME 


STATE 




OUTPUT 


iT 

XT ( 1 ) 


XT ( 2 ) 


XT ( 3) 


DELTC 


n.o 

o. loooooon 

01 

0.0 


0.0 


-0.75462280 

00 

0.01 

0.99640410 

00 

-0 . 1 673903D 

00 

-0.49775060-01 

-0.78731170 

00 

0.0? 

0.99119990 

00 

-0.32773581) 

no 

-0. 98343630- 

■01 

-0.R16825RD 

00 

0,03 

0 . 98446230 

00 

-0.48093550 

00 

-0.14558150 

00 

-0.84325030 

00 

0.04 

0.97626680 

00 

—0 .62691590 

00 

-0.19137800 

00 

-0.86667360 

00 

0.05 

0.96668910 

00 

-0.7656292D 

00 

-0.23563540 

00 

-0.88718710 

00 

0.06 

0.95580510 

00 

-0.89705180 

00 

-0.27826850 

00 

-0.90488470 

00 

0.07 

0 .94369090 

00 

-0.10211830 

01 

-0.31920330 

00 

-0.91986230 

00 

0.08 

0.93042170 

00 

-0. 11380420 

01 

-0.35837720 

00 

-0.93221780 

on 

0.09 

0.91607260 

00 

-0. 1247671D 

01 

-0.39573790 

00 

-0.9A205020 

00 

0.10 

0.90071780 

00 

-0.13501260 

01 

-0.43124270 

00 

-0.94945930 

00 

0.11 

0 * 88 44306D 

00 

-0. 14454850 

01 

-0.46485 84D 

00 

-0.95454750 

00 

0.12 

0.86728340 

00 

-0.15338380 

01 

-0.49656040 

on 

-0.95741470 

00 

0.13 

0.84934720 

00 

-0. 16152920 

01 

-0.52633220 

00 

-0.95816290 

00 

0.14 

0.83069190 

00 

-0.16899660 

01 

-0.55416490 

00 

-0.95689360 

00 

0.15 

0.8)138570 

00 

-0.17579910 

01 

-0.58005680 

00 

-0.95370790 

00 

0.16 

0.79 14954n 

00 

-0.18195110 

01 

-0.60401260 

00 

-0.94870640 

00 

0.17 

0.77108610 

00 

-0. 1874678D 

01 

-0.62604350 

00 

-0.94198870 

00 

0.18 

0.75022120 

00 

-0.19236540 

01 

-0.64616620 

00 

-0.93365360 

00 

0.19 

0.72896230 

00 

-0, 19666 1 OD 

01 

-0.66440250 

00 

-0.92379840 

00 

0.20 

0.70736890 

00 

-0.20037200 

01 

-0.68077930 

00 

-0.91251950 

00 

0.21 

0.68 549870 

00 

-0.20351700 

01 

-0.69532770 

00 

-O.R9991120 

00 

0.22 

0.66340740 

00 

-0.20611470 

01 

-0.70808310 

00 

-0.88606650 

00 

0,23 

0.64114880 

00 

-0.20818440. 

01 

-0.71908430 

00 

-0.871 0764D 

00 

0.24 

0.61877430 

00 

-0.20974580 

01 

-0.72837340 

00 

-0.85503010 

00 

0.25 

0.59633340 

00 

-0.21081871) 

01 

-0.73599580 

00 

-O.R3801440 

00 

0.26 

0.57387340 

00 

-0.21142340 

01 

-0.741 9991D 

00 

-0.82011440 

00 

0.27 

0.55143960 

00 

-0.21158030 

01 

-0.74643340 

00 

-0.80141250 

00 

0.28 

0.52907510 

00 

-0.21130960 

01 

-0.74935100 

00 

-0.78198900 

00 

0.29 

0.50682060 

00 

-0.21063190 

01 

-0.75080560 

00 

-0.761921R0 

00 

0.30 

0.48471510 

00 

-0.20956760 

01 

-0.75085250 

00 

-0.74128620 

00 

0.31 

0.46279520 

00 

-0.20813710 

01 

-0.74954R4D 

00 

-0.72015510 

00 

0.32 

0.44109540 

00 

-0.20 636060 

01 

-0.74695080 

00 

-0.69859870 

00 

0.33 

0.41964820 

00 

-0.20425820 

01 

-0.74311780 

00 

-0.67668480 

00 

0.34 

0.39848410 

00 

-0.20184970 

01 

-0.73810830 

on 

-0.65447840 

00 

0.35 

0.37763140 

00 

-0. 19915470 

01 

-0.73198130 

00 

-0.63204200 

00 

0.36 

0.35711640 

00 

-0. 19619260 

01 

-n. 72479610 

on 

-0.60943530 

00 

0.37 

0.33696360 

00 

-0. 19298220 

01 

-0.71661190 

00 

-0.58671560 

00 

0.3R 

0.31719550 

00 

-0.18954241) 

01 

-0. 70748750 

00 

-0.56393720 

00 

0.39 

0.29783270 

00 

-0. 18589130 

01 

-0.69748171) 

00 

-0.54115220 

00 


(c) Output — Concluded. 
Figure 3.- Concluded. 



The user’s main program— Some of the details of the main program are discussed briefly. The 
various matrices are first dimensioned and stated to be double precision. The problem will be solved 
using basically 3X3 matrices, but Z and PHI are 6 X 6 matrices so MAXRC is set to 36 (line 6). A 
double size dummy array is required in ETPHI, so DUMMY is dimensioned at 72, and KDUM is set 
to 72 (line 7). 

In line 8 the timing information is read in. TT is the initial time, DELT1 is the time increment 
used in the computation of P, DELT2 is the time increment, r j , desired in the printout of the 
transient and TFINAL is the final time for the transient. 

Lines 9 and 10 read data cards to fill a total of 7 matrices. 

Lines 14, 15, and 16 set up the appropriate constants for RICAT, specifying a print every step 
(line 14), the maximum number of steps to be taken by RICAT (line 15), and that both P and K 
should be printed (line 1 6). Lines 1 7 and 1 8 store the initial values of x 0 and P 0 in the running 
matrices, and lines 19 through 23 do the necessary computations to obtain P and K (called CK in 
program). Then F* and the transition matrix A2 (lines 28 through 32) are computed and printed. 
The transition matrix is labeled on the output (lines 29 through 31). Lines 33 through 39 page the 
output, print a heading for the transient response, and print the first point. Lines 40 through 47 
then increment the solution and the time, and print x(t) and u(t) (called XT and DELTC in the 
program). 


Example 4 — Sampled Data Ricatti Solution 

This example is provided to show the general use of the subroutine SAMPL. The theory of the 
example is given in the ASP manual, page 222, and very briefly in the dictionary description of 
SAMPL, page 24, in this manual. A listing of the main program is shown in figure 4(a). The data 
deck is shown in figure 4(b), and the output in figure 4(c). 

The main program is reasonably self-explanatory. The statement NCONT(2) = 4(line 13) 
indicates that SAMPL is to compute P for four successive time intervals and then stop. Both P and 
K (line 14) are to be printed at every step (line 12). 

As mentioned in the dictionary, K is the weighting matrix corresponding to the beginning of 
the interval, and P is the covariance matrix corresponding to the end of the interval. This is appar- 
ent in the output. For example, the first entry to SAMPL prints step number 0 and the K matrix, 
followed by step number 1 and the P matrix. On exit from SAMPL, P and K contain the data 
corresponding to Pj and Kj_! , which is the last interval. If printing is requested, the exit value of P 
and K will always be printed, and will be the last set of data. 
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c 


CHECK PROCEDURE FOR SAMPLE SEE PAGES 234 AND 244 OF ASP MANUAL 


0001 
non o 

DIMENSION NPH I ( 2 ) , NQ ( 2 ) ,NR ( 2 ) , NP ( 2 ) » NK ( 2 ) »NC0NT (3 ) »NH (2 ) 

nniiBi c Dccr ictnw duti- 2 ii u i i 

0003 
nnn a 

1 DUM ( 54 ) 

COMMON /MAX/MAXRC 

M AYDT-O 

— U 

0005 

0006 

— -Pf W A-n 

MDUM=54 

— U trtW 

0008 

0009 

LA L L KL) 1 1 \ L 

10 CALL READ ( 4, PHI , NPHI t H, NH , 0 ,NQ ,R ,NR ,R ,NR ) 

NP ( 1 ) =NPH I ( 1 ) 

lun / T \ uni i f / n t _ - ... 

— UU 1U 

0011 
0012 
nm -a 

IMP \ c ) =Nr n M £ ! 

CALL UNITY <P,NP) 
NC0NT ( 1)=1 

0014 

0015 
nn i a 

NC0NT ( 3) =4 

CALL SAMPL ( PHI , NPH I , H, NH , Q, NQ, R ,NR , P ,NP ,K ,NK ,NC0NT , DUM ,NDUM ) 

fAII FVTT .. .. 

— UU i. o 

0017 

0018 

UMUL — Er A -A- t — 

GO TO 10 
END 


(a) Main program. 


TEST 

PHI 

0 

PROGRAM FOR 
3 3 

1.0 

SAMPL CASE 1 FROM ASP MANUAL P234 AND P244 

0 

0 


0 

0 

0 


2 .0 

H 

2 

3 


0.0 

2.0 


0.0 

0.0 

0.0 


1.0 

0 

3 

3 


3.0 

1.0 


0.0 

1 .0 

1.0 


0.0 

0 

0 


1.0 

R 

2 

2 


1 .0 

1.0 



1.0 

2.0 




(b) Data. 

Figure 4.— Example 4. 



TEST PROGRAM FOR SAMPL CASE 1 FROM ASP MANUAL P234 AND P244 


VASP PROGRAM 


PHI 

0 j.0 

MATRIX 

3 ROWS 

_J nnnnnnm-v an 

3 COLUMNS 

no 


0.0 

0.0 


0.0 

0.0 


0.0 

2.0000000D 00 


H 

0.0 

— M 

MATRIX 

2 ROWS 

2.00000000 00 
-GU-0 

3 COLUMNS 

0.0 

l.OOOOOOOD 00 



0 MATRIX 3 ROWS 3 COLUMNS 

3.000 0 0000 00 l.OOOOOOOD 00 -OrO 

1 .00000000 00 l.OOOOOOOD 00 0.0 

0.0 0.0 l.OOOOOOOD 00 


R MATRIX 2 ROWS 2 COLUMNS 

l.OOOOOOOD 00 l.OOOOOOOD 00 

— l.OOOOOOOD 00 2 . ■ 0 0 000000 -00 

STEP NUMBER= 0 IN SAMPL 


KIT) MATRIX 3 ROWS 2 COLUMNS 

4.2857143D-01 -1 .4285714D-01 

— 0^0 OrO 

-1.4285714D-01 7. 1428S71D-01 

STFP Nl IMRFR = 1 TN RAMPI 

PIT) MATRIX 3 ROWS 3 COLUMNS 

— 3 . 14285710 OQ l.OOOOOOOD 00 2 . 8571 4 29D-01 

l.OOOOOOOD 00 l.OOOOOOOD 00 0.0 

2 .85714290-01 0.0 3.5714286D 00 


(c) Output. 
Figure 4.- Continued. 



STEP NUMB ER= 1 IN SAMPL 

K ( T ) MATRIX 3 ROWS 2 CPI UMNS 

4. 1489362 D— 01 -7 .4468085D-02 

0.0 0.0 

t2 . 65957 4 50-01 1 . 3297872D 00 

STEP NUMBER* 2 IN SAMPL 

PIT) MATRIX 3 ROWS 3 COLUMNS 

3.1702128D 00 l.OOOOOOOD 00 5.31914890-01 

l.OOOOOOOD 00 l.OOOOOOOD 00 0.0 

5.31914890-01 £U-Q 5 . 78723 4 QO 00 

STEP NUMBER* 2 IN SAMPL 


KIT) MATRIX 3 ROWS 2 COLUMNS 

4.10 54403D-0 1 -5.2720135D-02 

-■■■ 0 .' Q— — — ■ 0.0 — — 

-3.0510376D— 01 1.52551880 00 

STEP NUMBER* 3 IN S A-MPL 

PIT) MATRIX 3 ROWS 3 COLUMNS 

3.1780119D 00 liO O OOOOOD 00 6*1020752D 01 

l.OOOOOOOD 00 l.OOOOOOOD 00 0.0 

6. 1020752D-01 0.0 6.4918676D 00 


STEP NUMBER* 3 IN SAMPL 

KIT.) MATRIX 3 ROWS 2 COLUMNS 

4.0964801D-01 -4. 8240037D-02 

0.0 0.0 

-3 . 13167930—01 1 . 5658397D 00 

STEP NUMBER* 4 IN SAMPL 

PIT) MATRIX 3 ROWS 3 COLUMNS 

3.1 807040D 00 l.OOOOOOOD 00 6.2633587D-01 

l . OOOOOOOD 0 0 U Q Q QOOOOD 00 -Q^Q 

6.2633587D-01 0.0 6.6370228D 00 

(c) Output — Concluded. 
Figure 4.— Concluded. 



Example 5 — Matrix Decomposition 


This example is a test program to check the operation of DECGEN. It first generates a matrix 
R to be decomposed, then proceeds with the decomposition, and checks the result, printing all of the 
associated matrices. The general procedure is to input a diagonal matrix ZL and transform it into the 
matrix R to be decomposed. Figure 5(a) is a listing of the main program; figure 5(b) is a listing of 
the subroutine ORTH; figure 5(c) is the data deck; and figures 5(d) through 5(f) are the output. 

In the main program, all matrices are dimensioned 100, although the actual matrix size used is 
2X2 and 4X4. Accordingly, MAXRC is set to 100. The dummy matrix is dimensioned 700, 
since DECGEN requires that much. The input matrices ^re read at line 8. 

Subroutine ORTH, called at line 9, produces a n X n orthogonal matrix, using the original T 
matrix, and places the results back in T. The procedure is as follows. 

First, generate an elementary rotation matrix Ejj. This is a unity matrix, with elements eii and 
ejj replaced by cos tjj and elements ejj = —ejj = sin tjj. 

Then, 

T = n Ejj 

Lines 10 through 17 set up indices for referring to the seven dummy matrices. The input 
matrix, ZL, is then transformed by the matrix T, so that 

ZL, = T*ZL*T' 

Note that ORTH leaves T' in DUM3. Also, if the T at input was the null matrix, the rotation will 
be the identity matrix, so that R = ZL. Lines 19 through 27 then juxtapose either the matrix EXR 
or the matrix EXC, using JUXTR or JUXTC, depending on the compatibility of the dimensions. If 
both sets of dimensions are incompatible, no juxtaposition is done. In any case, the result of this 
operation is placed in R. The decomposition routine is called next. If the original ZL matrix had 
zero in element (2,1) and no juxtaposition was done, then R is assumed symmetric, and the 
DECSYM entry is used. If ZL was not symmetric, the program will produce errors. Otherwise, 
the DECGEN entry is used (lines 29 to 31). Finally, the resulting matrices H and G are tested 
using 


R, = HG 
RE = R - R, 


and all resulting matrices are printed. 

In figure 5(c), blank lines represent blank cards. In the data cards for case 4 the header card 
for EXR has no dimension information and no associated data cards. This indicates that the 
matrix EXR is to be left unchanged, and that no data cards are to be read for EXR. In case 7, 

EXR is again left unchanged. A blank data card follows the EXC header card. 

The output (figs. 5(d) through 5(0) contains the results of decomposing three different matrices. 
Figure 5(d), case 1 , is a 2 X 2 rank 1 matrix; figure 5(e), case 4, is a 2 X 3 rank 2 matrix; and 
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-P> 


0001 


MAIM PROGRAM TO CHECK DECOM ET AL 

DIMENSION NZL ( 2) ,NT( 2) , NEXRI 2) ,NR(2) ,NG(2 ) ,NH(2) ,ND(2) ,NR1(2) , 


0002 
r\ r\r\~% 

1 

1 

NRC(2) fNCXCm 

DOUBLE PRECISION ZL ( 100 ) , T ( 100 ) , EXR ( 1 00 ) , R ( 100 ) , G ( 100 ) , H ( 100 ) , 
DUM(700) ,R1( 100) ,RE( 100) ,EXC( 100) 

r' oil u ah /u t w j u i w n r ... 

UUU J 

0004 

0005 

nnrv l 


LGMnUN /“AX/MAXRL 
MAXRC = 100 

0007 

0008 
OOOQ 

20 

RUUM 1 * t OO 

CALL RDTITL 

CALL READ ( 4 , Z L , NZ L , T , NT t EX R, NE XR , E XC ,NEXC , T»NT ) 

r a i i no th / t mt riim i/nitM ^ _ _ _ 

0010 

0011 

0012- 


l/ALU tmitl r i fINf fUUr1fT\lJUrn I 

M»NT( 1) 

M2=M*M+1 

M C-M*IU1 ...... . . 

0013 

0014 
noi^ 


M3=M2+MS 

M4=M3+MS 

0016 

0017 

Wllfl.. 


M6=M5+MS 

M7=M6+MS 

nilM(M7 \ m 7 1 19 V _ 

0019 

0020 
on? l 


CALL MULT ( T , NT , ZL , NZ L, DUM ( M3 ) , ND ) 

CALL MULT ( DUM ( M3 ) , ND, DUM ( M2 ) , ND , ZL ,NZL ) 

TC /M71 / II 1IC KlCV/TMUOn TH AA __ _ _ 

Ulft I 
0022 
0023 

Amy 

<■* A 

ir INiL 1 1 ; •Nt^lNtAH IMuU IU 30 

CALL JUXTC < ZL»NZL, EXC, NEXC , R»NR) 
GO TO 50 

0024 

0025 

0026 
r\ r\ n -r 

30 

t A 

IF (NZL(Z) « NE • NE aR ( 2 ) ) GO TO *+0 
CALL JUXTR (ZL,NZL,EXR,NEXR,R,NR) 

GO TO 50 

A i 1 i n r\i i » *r r* < • i I i I r-» 1 ■ n i . . . .. ...... ...... ......... .... 

■0027 

0028 

0029 

-00-36 

r 

CALL EOUAT E ( ZL , NZ L , R ,NR ) 

IF (DUM(M7) .NE.O.DO) GO TO 50 

CALL DECSYM (^rNRf,G , NG ,H , NH , DUM , KDUM ) 

&e -TO 76 - 


(a) Main program. 
Figure 5.- Example 5. 



0031 50 CALL DECGEN ( R t NR ,G,NG,H f NH,DIJM,KDUM) 

0032 70 CALL MULT ( H*NH,G»NG,R 1,NR1 ) 

QQH CALL SU B T (R -T NR T RlyNM - T^ E yNR E - ) 

0034 CALL PRNT (R,NR,'R i,l) 

0035 CALL PRNT (R1,NR1,'R1 M) 

CALL ft RMT ( RE T W5 , » R E RR 1 , 1 ) 

0037 CALL PRNT (H ,NH ,'H *,1) 

0038 CALL PRNT (G ,NG »'G M) 




0040 

NU \ 1 ) = Jr 

ND ( 2 ) = 1 


0041 

0042 

-0G43 

CALL PRNT ( DUM ( M6 ) * ND » ' RANK ' » 1 ) 

GO TO 20 

frNB 



(a) Main program - Concluded. 
Figure 2.- Continued. 


0001 

0002 
— 000-3 

0004 

0005 
— 0006- 

0007 

0008 
— 0009 

0010 

0011 
— 00 42 

0013 

0014 
0045- 

0016 

0017 
0048- 

0019 

0020 
0044 

0022 

0023 
0024- 

0025 

0026 
0027 

0028 

0029 

—0030 

0031 


SUBROUTINE ORTH ( T , NT , DUM , KDUM ) 

DOUBLE PRECISION T ( 1 ) , DUM ( 1 ) ,CTH, STH 

DIMENSION NT ( 2 ) 

LDM2=2*NT( 1)**2+1 
10 LDM1=NT(1)**2+1 

N = MT( 1 - ) 

CALL UNITY ( DUM ( LDM1 ) , NT ) 

NM=NT ( 1 ) — 1 

DO 20 J-1 ,M M 

JP=J+1 

DO 20 I=JP,N 

CALL UNITY (BU M, - NT T 

1 1 =N* ( I -1 ) + 1 
JJ=N*( J-l )+J 

1 - 1 KJ 

JI=N*( J-l)+I 
CTH=DCOS(T( I J) ) 

ST II- DS IN ( T 

DUM (II) =CTH 
DUM ( J J ) =CTH 

DU M ( I J ) -ST H 

DUM ( J I )=-STH 

CALL MULT ( DUM ( LDM 1 ) , NT, DUM , NT, DUM ( LDM2 ) , NT ) 

-OAtt " E OU AT T: ( DUM( L D M 2 ) , NTrOU M ( LDM 1 ) , NT ) 

20 CONTINUE 

CALL TRANP (DUM(LDMl) ,NT,T,NT) 

CAL fc -M UL - T -- (Tf f H ^ DU M ILDMl ) r NT ( LP tt g ) ,NTU — 

CALL PRNT ( T , NT , 4HT ,1) 

CALL PRNT ( DUM ( LDM2 ),NT,4HT*T',1) 

RETURN — 

END 


(b) Subroutine ORTH. 
Figure 5.— Continued. 



TEST 

PROGRAM 

FOR DECGEN AMO DECOM 

CASE 1 

2X2 RANK 1 

ZL 

2 

2 



1.0 

1.0 




2.0 

2.0 




T 

2 

2 




EXR 

1 

1 






EXC 

1 

1 






TEST 

PROGRAM 

FOR 

DECGEN 

AMD 

DECOM CASE 

4 

2X3 RANK2 

71 

2 

? 






1. 

2.0 







T 

2 

2 





- 


.7 







EXR 








EXC 

2 

1 






2. 








3. 








TEST 

PROGRAM 

FOR 

DECGEN 

AND 

DECOM CASE 

7 

I LL-COND 4X4 RANK 3 

ZL 

4 

4 






1. 









2. 




l.D-6 



T 

4 

4 







.2 


.3 


.4 






.5 


.6 








.7 



EXR 








EXC 

1 

/ 







(c) Data. 

Figure 5.— Continued. 



2X2 RANK1 


V ASP PROGRAM 


oo TEST PROGRAM FOR DEC GEN AND DECOM CASE 1 


ZL MATRIX 2 ROWS 2 COLUMNS 

lr S OOOO&OO-OO i*OOOOOOOD 00 

2 • OOOOOOOD 00 2 . OOOOOOOD 00 

T MATR - I -X 2— ROWS 2-GO LUM W S 

0.0 0.0 

0.0 0.0 


EXR MATRIX 1 ROWS 1 COLUMNS 

0.0 


EXC MATRIX 1 ROWS 1 COLUMNS 

0.0 


T MATRIX 2 ROWS 2 COLUMNS 

1. OOOOOOOD 00 0.0 

XUD 1. OOOOOOOD 00 

T*T ' MATRIX 2 ROWS 2 COLUMNS 

1 . OOOOOOOD 00 0^0 

0.0 1. OOOOOOOD 00 

R HATR-LX 2_RQWS 2~ -COLUMNS 

1. OOOOOOOD 00 1. OOOOOOOD 00 

2. OOOOOOOD 00 2. OOOOOOOD 00 


R1 MATRIX 2 ROWS 2 COLUMNS 

1. OOOOOOOD 00 1. OOOOOOOD 00 

-2-.-6 00000C H3 — DO -- ~~2~.~ 00 00000 B ~ (TO 

RERR MATRIX 2 ROWS 2 COLUMNS 

4 1 1 6 333 6 3D 10 4 .1633363D 16 

4.4408921D-16 4 . 4408 92 1 D-16 


(d) Case 1 . 

Figure 5.— Continued. 



H-..- MATRIX- 2 ROWS- - — 1--GML44M-NS 

1.41421 36D 00 
2 .8284271 D 00 


G MATRIX 1 ROWS 2 COLUMNS 

7.0710678D-01 7.07 10678D-01 


RANK MATRIX 1 ROWS 1 COLUMNS 


1 . OOOOOOOD 00 


(d) Case 1 — Concluded. 
Figure 5. - Continued. 


TEST PROGRAM FOR DECGEM AND DECOM CASE A 


2X kaNK2 


V ASP PROGRAM 


ZL MATRIX 2 ROWS 2 COLUMNS 

l . OOOOOQOD 00 0^0 

0.0 2 • OOOOOOOD 00 

1 MATRIX 2 ROWS 2— C OLU M NS 

7. 99569850-01 6.0057311D-01 

-6.00573110-01 7.9956985D-01 


EXR MATRIX 1 ROWS 1 COLUMNS 

0.0 


EXC MATRIX 2 ROWS 1 COLUMNS 

2. OOOOOOOD 00 

3.00000000 00 

T MATRIX 2 ROWS 2 COLUMNS 

8 . 2501188D-01 -5 . 65115390 - 01 

5 • 6511 539D-01 8 .2501 188D-0 1 

I&XI MATRIX 2 -ROWS 2 COLUMNS 

1. OOOOOOOD 00 0.0 

0.0 1. OOOOOOOD 00 


R MATRIX 2 ROWS 3 COLUMNS 

1.3193554D 00 -4.66226910-01 2.00000000 00 

-4 . 6 6 226910-01 1 . 6806 44 6D 00 3 . 000000(30 00 


R 1 MATRIX 2 ROWS 3 COLUMNS 

-4,6622691 D-0 1 1.6806446D 00 3. OOOOOOOD ''O 


(e ) Case 4. 

Figure 5.— Continued. 





RE££t MATRIX 2 RO WS 3 CQU JMNS - 

2 • 22044 60D-16 -2 .49800180-16 4.4408921D-16 

-9.7144515D-17 2 . 2204460D-16 4 . 440892 1 D-16 


H MATRIX 2 ROWS 2 COLUMNS 

2 .04935730 00 1.3259717D 00 

OiO 3 * 4 70 1 4 90D 00 

G MATRIX 2 ROWS 3 COLUMNS 

7 w 307 190 6 D - 0 1 - 5* 4 0 8 5 9 65 0 - 01 4.-1655791-0-6-1- 

-1.3435357D-01 4.84314830-01 8 . 6451620D-01 


RANK MATRIX 1 ROWS 1 COLUMNS 

2 . OOOOOOOD 00 


O' 


(e) Case 4 — Concluded. 
Figure 5.- Continued. 



TEST PROGRAM FOR DEC GEN AND DECOM CASE 7 ILL-COND 4X4- RANK 3 


VASP PROGRAM 


ZL MATRIX 

4 ROWS 

4 COLUMNS 


i -nnnnnnnn nn 

o*o 


QO- 

0.0 

2.00000000 00 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Q.Q 

-CLrO 

-0-.0 

1 » 00000000-06 


T MATRIX 4 ROWS 4 COLUMNS 


8.7748949D-01 

-3.9329146D-01 
-2.6426304D-01 • 

-7 /1CC7^0n-AO 

1 .54979 ITD-Ol ■ 
7.6892918D-01 
-5.4846577D-01 

-7 QOZOA7A r\ — f\ 1 

2*40715210-01 
2. 066722 3D-01 
6.9767211D-01 
a Ann-m 

3.75650570-01 

4. 59735070-01 
3.77629410-01 

f lyy t 37U*Uc — 

EXR MATRIX 

0.0 

“c «t5700UJU U \J 1 

1 ROWS 

•0 07c0 lOUU"U 1 

1 COLUMNS 

f r?U"Ul 


EXC MATRIX 1 ROWS 1 COLUMNS 

^0 


T MATRIX 4 ROWS 4 COLUMNS 

8 .8 94 8 5 4 8 0 - 01 1.3 8 9 56 7 2 0 - 01 2*20630900 - 01 

-6.6804271D-02 8.9817561D-01 -1 .3776143D-0 1 

1.1688667D-01 1 . 2733255D-02 9. 4444935D-01 

4. 36 8 031 6 D - Q1 4 . 1 6 90 464 D - 01 1* 9 17 42 2 -6 0 -0-1 


- 3.705 9 5760 - 01 
-4. 1211595D-01 
-3.0690520D-01 
7.7371082D - 01 


T*T* MATRIX 4 ROWS 4 COLUMNS 

- 1 . Q O 000000-0 0 2 w - 7 - 7 5 5 57 -6 0— 1 - 7 1 *3 8 7 - 7 - 7 8 00 - 1 - 7 - 

2 .7 75 5576 D- 17 l.OOOOOOOD 00 4. 1633363D-17 

-1.38777 88D-17 4. 1633363D-17 l.OOOOOOOD 00 

-5 . 551 1 151 D- 17 4 . 1 6 333 6 3D — 17 0^9 


1 5 « • 55 1 1 1 5 t ~ B — 1 - 7 
-4.1633363D-17 
0.0 

1 . 0000000 0- 00 


(0 Case 7. 

Figure 5.— Continued. 



R MATRIX 4 ROWS 4 COLUMNS 


- 8 . 2 9 6 9 5 770 - 01 3^09 032 3 4 0 - 0 1 

-3. 09032340-01 1.6179018D 00 

1.0042376'>-01 1.5064996D-02 


1. 00423300 - 01 
1. 50649960-02 
1 .398 6860D-02 


2.72 640 1 8D-0 1~ 

7.19726510-01 
6. 1673338D-02 


R1 MATRIX 4 ROWS 4 COLUMNS 

0.2969 5 770 - 01 3 i090323 ' 4D - 0 1 1.0042336D - 01 

-3.09032340-01 1.6179018D 00 1 . 506499 60-02 

1.0042336D-01 1.50649960-02 1.39868010-02 

- - 2.726401 8D-01 7 ri9- 7 2651D - -01 67 167362 8D-02 


2. -7 2 - 6401 8 D- 0 1 
7 • 1972651 D-01 
6.16736280-02 
5. 38 4 15120-01 - 


RERR MATRIX 4 ROWS 4 COLUMNS 

6.9388939D-17 — " 8 . 32 66 7270 - 17 A » 1 6 333630 —37- 

-5.55111510-17 2 • 2204460D-16 -1.73472350-18 

4 . 1 633363D-17 0.0 5 . 8972036D-08 

— 1 . 38777 88D-17 1.5265567 0- 1 -6 2-. 90 4 7 379D - 07 


1 i 38777000 - 17 
8.32667270-17 
-2.90473790-07 
— 1. 4 30 7633D -06- 


H MATRIX 4 ROWS 2 COLUMNS 

8 .77 8 7704D - 01 2. 4 20 56- 120 - 01 

0.0 1.27196770 00 

1.17671260-01 1 . 1843851D-02 

4 . 6 71 65 390 - 01 5 .6 5 837100 - 01 


G MATRIX 2 ROWS 4 COLUMNS 

8 . 770770 6D-01 &t 0 1 . 1767126D -0i 4. 671 6539D - 01 - 

-2. 42956120-01 1.2719677D 00 1.18438510-02 5 . 6583710D-0 1 

— — RANK — M ATRI -X-- L-RGW5 COLUMNS 

2. OOOOOOOD 00 


O' 

CO 


(f) Case 7 — Concluded. 
Figure 5. - Concluded. 



finally, figure 5(f), case 7, is a 4 X 4 matrix of rank 3, with one very small eigenvalue equal to 1 O' 6 . 
The error matrices of the first two decompositions are extremely small, but that from the third one 
has errors of the order of 10‘ 6 . These are caused by the built-in pivot rejection device, which 
rejects all pivots smaller than 2X1 0“ 5 times the largest of the diagonal elements (see DECOM, p. 85 
and PSEU, p. 70). This last matrix, case 7, was also tried with an eigenvalue of 10' 3 , and the errors 
were then on the order of 1 0 _1 6 . 


Example 6 — Use of the Pseudoinverse Routine 

This program is designed to check the operation of PSEUDO. The procedure is as follows: 

First the input matrix A is read; then B = A# is computed. The accuracy of the pseudoinverse 
is then checked by the first two Moore-Penrose axioms 

BAB - B = A e 

ABA - A= B e 

All the various matrices are printed. 

Figure 6(a) is the program listing and figure 6(b), the output. Three cases are presented; the 
first two are the examples presented in the ASP manual; and the third one contains several zeros. 

The first matrix printed for each case is the input matrix and each has a different label. The other 
titles are abbreviations chosen to fit the allotted four character space as follows: 

APSE -»■ A # 

AASA -* AA # A 

AERR -> Ap or B £ 

ASAA -»■ A# A A# 

It can be noted that the size of the numbers in the AERR matrices is 1 0‘‘ 6 , which is very good. 

Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., 94035, July 12, 1971 
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0001 


0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 

- 0011 

0012 

0013 

0016 

0015 

0016 

- 004-7 

0018 

0019 

0020 

0021 

0022 

-0023 

0024 


DIMENSION A (50) , B ( 50 ) ,W(350) ,NA( 2), NR (2) f A1 (50) ,A2 ( 50) *NA1 ( 2) 
1 NA2 ( 2 ) » L AB ( 2 ) 

DQU OL E - PRECISION A t B t W t Al y A 2 

COMMON /MAX/MAXRC 


5 CALL RDTITL 

CALL READ ( 1 , A, NA » A , NA» A , NA , A, NA , A , NA ) 

CALL PSEUDO( A,NA,B,NB,W,NW) 

CALL PR NT (B,NB, 'APSE',1) 

- C A L L . MULT( A» N Af-&^Na^A4^NA],l 
CALL MULT( A 1 ,NA 1 , A, NA, A2 , NA2 ) 

CALL SCALE( A , NA, A 1 , NA 1 , -1 . DO ) 

■ CA LL ADD ( A 1 ,- N A l -r- A2 , NA2 , A 1 , NA 14-- - - — 

CALL PRNT ( A2,NA2, • AASA' , 1 ) 

CALL PRNT (A1,NA1, 'AERR* ,1) 

r ai i Mi ll T ( R . NR . A . N A . A 1 . N A 1 ) -- - - 

CALL MULT( A 1 ,NA 1 , B,NB, A2 , NA2 ) 

CALL SCALE(B,NB,A1,NA1,-1.D0) 

CALL ADD(Al ,N Al r A2 , NA2 ,ArT:-rNArl»- 

CAlL PRNT <A2,NA2» 'ASAA',1) 

CALL PRNT ( A1 » NA 1 , ' AERR ' , 1 ) 

. ..oo-m 5 - - - - 

END 


(a) Main program to check PSEUDO. 


Os 

L /1 


Figure 6.- Example 6. 



On 

On 


PSEUDO TFST PROGRAM CASE 1 


FROM ASP MANUAL PAGE 146 


VASP PROGRAM 


B MATRIX 3 ROWS 4 COLUMNS 

4 . 00000000 00 1 . OOOOOOOD 00 3.0 0 00000 0 00 £- .00000000 0 0 

-2.0000000D 00 5. OOOOOOOD 00 -1. OOOOOOOD 00 -3. OOOOOOOD 00 

2. OOOOOOOD 00 1 .3000000D 01 -9. OOOOOOOD 00 -5. OOOOOOOD 00 


APSE MATRIX 4 ROWS 3 COLUMNS 

9 . 5029697D-02 -5.6580181D-02 2 .031R850D-02 

-6 .736480 20-02 3 . 1 884964D-02 -3.90747UD-02 

5.0640825D-02 -3 .673022 8D-02 -8 .9090341 D-03 


AASA MATRIX 3 ROWS 4 COLUMNS 


4. OOOOOOOD 00 -1. OOOOOOOD 00 -3. OOOOOOOD 00 2. OOOOOOOD 00 

* 2 . 00000000 -- aa... 5 ^ooqqqqoo..oq. . .*- 1.00000000 oa - .= 4.00000000 oo 

2. OOOOOOOD 00 1 • 3000000D 01 -9. OOOOOOOD 00 -5. OOOOOOOD 00 


AERR MATRIX 3 ROWS 4 CQLLIMNS- 

-6. 6613381 D- 16 6 . 661 338 ID- 16 2 . 2204460D- 1 6 -6.6613381D-16 

1.3322676D-15 0.0 -8 . 88 17842D-16 6.6613381D-16 

~2-^2 Q446QQ=-- l-5 34996402 90=1 5 -^4^240847 50=4-5- - -6.661 33410-1-6 


ASAA MATRIX 

9.5029697D-02 
-3.0790872D-02 
-6 . 7364802D-02 
. .5436 40 8 2 5 0=02 - 


4 ROWS 

-5.65-801810-02 
3.31 35355D-02 
3.1 8849640-02 
*4*42 302240=02- 


3 COLUMNS 

2.03188500-02 
3. 78243200-02 
-3. 90747110-02 
.. ==84*090341 0*03 


AERR MATRIX 4 ROWS 3 COLUMNS 


"1.38777880-17 
2 . 341 8767D-17 
0.0 

-1.82 145940- 1 7 


1.21430640-17 
-1 .04083410-17 
-8. 6736174 D- 19 
7.8062 5360—18 


3.46944700-18 
1 . 6479873D- 1 7 
- 1 . 4745150D-17 
—4. 33680870- 1 8 


(b) Output. 


Figure 6.— Continued. 



PSEUDO TFST PROGRAM 


MATRIX 


CASE 2 
4 ROWS 


l.OOOOOOOD 00 
-2 . OOOOOOOD 00 
=£,,00000000 -00 


— i-^ooooGOOo-oo-- 
2 .5000000 D 01 
-fi. OOOOOOOD 00 
6* 000000 00 00 


APSE MATRIX 


ROWS 


-1 ♦ 1026095D-03 
-4.6911023D-02 ■ 

-^32 83103 0^02 - - 

AASA MATRIX 


2.9737044D-02 
-7.551 2045D-03 
. . <^75642360-03 

< 4 ROWS 


l.OOOOOOOD 00 2 • 5000000D 01 
■2. OOOOOOOD 00 -8. OOOOOOOD 00 
4*00000000 00 - 6*00000000 00 


AERR MATRIX 


4 ROWS 


■ 4 .44 0892 10 ^ 1 6 -- 0*0 

-4.4408921 D-l 6 -3 . 5527137D-15 
6.661 3381 D-l 6 6 .66 1 338 1 D- 1 6 
^22064606-46 -2*22 0 44 600 -16 


ASAA MATRIX 


ROWS 


-1 . 1026095D-03 
-4. 69 1 10 23D-02 - 

■6*32331030-02- 

AERR MATRIX 


2.9737044D-02 
-7.5512045D-03 
6*4 5-642 360-03 

C 4 ROWS 


0.0 0.0 

8 . 6736174D-1 8 1 .7347235D-18 

.1*3877.7880-17 .-2.60208520-18 


a\ 

-0 


VASP PROGRAM 




PSEUDO TEST PROGRAM CASE 3 


VASP PROGRAM 


MATRIX 


4 ROWS 


4W 0- 


-CWO- 


0.0 

3. 5000000 D-02 
-OO - 


-3.9 100000D 00 

0.0 

- 3.4-0000000-01 


2 COLUMNS 


APSE MATRIX 2 ROWS 4 COLUMNS 

TWO -3 . 1331 4 71 D - 02 5 . 501293 0 0 - 03 W9 5 1- 8 G 81- 0 -01- - 

0.0 -2. 557541 7D-01 2 .8046074D-04 3 . 879891 6D-06 


MSA.—- MAI RJLX. A -ROWS- 2 COL UMNS-- 

0.0 0.0 

-1.0889456D-17 -3.9100000D 00 

W.5OQ000QQt02 CUO 

-2 • 5300000D 00 3 . 1000000D-01 


AER.R MAIR.J-X 4- ROWS - 

0.0 0.0 

1.0889456D-17 4 . 440892 ID-16 

0^0 0^0 

2. 2204460D-16 0.0 


-ASAA MATR I-X ^ ROWS- A COLUMNS 

0.0 -3. 1331471D-02 5. 50129300-03 -3.95180810-01 

0.0 -2.5575417D-01 2 . 8046074D-04 3.87989160-06 


AERR MATRIX 2 ROWS 4 COLUMNS 

0.0 0.0 0.0 0.0 

-O^ 2 . 7755576D-17 5 ^ 4 2101090 - 20 1. - 101354 6 0 - 1 5 


(b) Output — Concluded. 
Figure 6.— Concluded. 



APPENDIX A 


DESCRIPTION OF INTERNAL SUBROUTINES 


25. READ1 
DESCRIPTION 

This subroutine reads a single matrix from cards, without a header card. It is called by READ, 
after the latter has read the header card. The dimensions of the matrix to be read are in array NZ. 
If this is zero, no array will be read. In any event, the routine then prints either the array just read, 
using NZ for dimensions, or, if NZ = 0, the array already stored, using NA for dimensions. 

The subroutine reads the data from cards, each row of the matrix starting on a new card, using 
format (8F 10.2). If the card data is in exponential form, it must use a D exponent. 

USAGE 

CALL READ 1 ( A,NA,NZ,NAM) 


Input Arguments 
Matrix: 

Dimension array: 
Constant: 


Output Arguments 
Matrix : 

Dimension array: 


A (if NZ = 0) 

NA,NZ 

NAM, containing a four-character (or less) name for the matrix, 
which will be used by PRNT 


A (if NZ =£ 0) 
NA 


26. ASPERR 
DESCRIPTION 

This is an installation dependent subroutine. It is called by the various subroutines when they 
detect an error. It is intended to provide an error walkback, so that the programmer can determine 
which call of a given subroutine is in error. It also counts the number of errors and calls EXIT 
after ten entries into ASPERR. 
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USAGE 


CALL ASPERR 

It has no arguments. The user may, if he wishes, call this program to help him track down errors. 

Subroutine ASPERR calls in turn a system program which provides the actual walkback. In 
Ames OS this system routine is called ERRTRA, while in Ames TSS, it is called TRACE. The 
calling statement should be changed to match the user’s operating system, or else deleted altogether. 


27. BLKDATA 

/ 

DESCRIPTION 

This is an installation dependent subroutine. It loads certain common areas used by VASP 
with appropriate constants as follows: 

1 . COMMON/F0RM/NEPR, FMT 1(6), FMT2(6) 

These three variables control the printing procedure, and are set to 7, (1P7D16.7), and 
(3x,lP7D16.7), respectively. They assume a line length of at least 1 15 characters. 

2. COMMON/ LINES/NLP, LIN, TITLE(23) 

NLP controls the number of lines per page, and is set at 45 to agree with the NASA— Ames 
system. It should be changed to match each installation. 

LIN is a counter which keeps track of the number of lines printed on each page. It is 
incremented and used only in LNCNT. 

TITLE contains 72 blank characters, which can be loaded as desired by use of RDTITL, plus 
20 more characters containing “VASP PROGRAM.” Subroutines LNCNT prints TITLE at 
the head of every page. 

3. COMMON/M AX/M AXRC 

MAXRC is used by most subroutines to check the reasonableness of the matrix dimensions. 
The user should set MAXRC to match the storage available for each matrix. It is preset to 
6400. 


28. PSEU 
SUMMARY 

PSEU is a FORTRAN routine to find thje Moore-Penrose generalized inverse of a non-negative 
definite double-precision matrix. It has a separate entry PSEUP for input of a matrix that is 
already symmetric. A symmetric matrix is always used for the actual diagonalization process. This 
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process is done in a self-contained subroutine, ANDRA. The routine “never fails, since it includes 
the singular case. However, it may fail to give the correct rank. To control this, an option to do side 
calculations is available. After the first pivots have been found, if the rank is not maximum, the 
result of each pivot step is used in two axiomatic expressions (subroutine BDNRM). This side cal- 
culation yields a measure of the worth of the pseudoinverse obtained so far. This result is multi- 
plied by a parameter factor raised to the power of the current rank (nonlinear penalty function). 

The routine can backtrack from the first bad step and stop with the previous rank. It has an option 
to do the minimum calculations for getting a rank only. The generalized inverse is useful for least- 
squares solutions of Ax = b; it works when A is singular. This method is best suited to symmetric 


matrices. The routine has suitable error exits. 
USAGE 

CALL PSEU(A,B,C,EE,DEP,IP,D) 
or 

CALL PSEUP(A,B,C,EE,DEP,ID,D) 
Note: PSEUDO usesPSEU entry. 
Input Arguments 

Matrix: A 

Control arrays: DEP 

DEP1 

DEP2 


Description 

The array, to be inverted, left intact, must be 
symmetric if PSEUP call is used. Non-negative 
definite, or nearly so. 

Values DEP1, DEP2, DEP3 

Default: If zero, user gets 2.D— 6 used instead. 
This number is multiplied times the largest mag- 
nitude on the diagonal of B at start. If any 
trial pivots are found less than this, they are 
avoided as zero. 

Default: If zero, user gets 1 .DO used instead. 
Needed only if iteration. The routine computes 
two numbers, p, q, which would be zero if the 
first two Moore axioms were satisfied. This num- 
ber is raised to the power of the number of 
pivots found as a factor to use to make the 
product with the sum of p and q larger. Mak- 
ing this product larger tends to make the routine 
reject the current pivot. Values between 1 and 
2 work for ordinary purposes. 

Note: PSEUDO uses default values of DEP1 
and DEP2. 
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DEP3 


Output Arguments 
Matrix: 

Matrix: 

Matrix: 


This is for output only. It holds the last pivot 
actually accepted. This gives the user or calling 
routine an estimate of the size of pivots found, 
in case effective rank is not that desired, operat- 
ing with given value of DEP1. If iterating, this 
may be the last pivot rejected. 

IP Parameter array of integers IP1, IP2, IP3, IP4. 

IP1 If zero, do not iterate with side calculations. 

If 1 , iterate. 

Note: Other values should not be used, since 
DECOM employs peculiar values. 

IP2 If zero, do all calculations, otherwise do rank 

only. 

Note: Setting this to zero for each call is very 
useful in avoiding confusion between ranks 
determined from different calls. Used also to 
output the effective rank. PSEUDO sets IP1 
and IP2 to zero. 

IP3 The row size of the matrix input. 

IP4 The column size of the matrix input. 

Note: IP4 need not be specified for PSEUP 
entry. 


B Holds the pseudoinverse output. (In rank only 

case, holds a diagonal matrix with 0’s and 1 ’s 
corresponding to pivots accepted or rejected.) 

C In nonsingular case, holds the matrix T of the 

diagonalization case. In singular case, holds that 
certain matrix U described in ASP manual. 

EE Holds the pseudoinverse of the original B. 

Note: A and B are the same size. The other 
matrices are square, of the size of C, which is 
determined by the smaller dimension of A. D is 
either five times the size of C, if iterating, or the 
same size as C. 
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Matrix: 


D 


In the nonsingular case, D holds a copy of the 
B formed from A. (It equals A for a PSEUP 
entry.) In the singular case, it holds a pseudo- 
inverse for a “B” permuted so that independent 
variables are all moved to the left-most positions. 
Note: D has possibly four other matrices. Let 
these be Dl, D2, D3, and D4, in order. They are 
used only if iterating (Dl also used by DECOM). 
Dl, D2 hold old results. D3, D4 holds intermedi- 
ate values when doing the side calculations. 
PSEUDO does not provide for Dl through D4. 


Notes on Usage 
Symmetry 

This method is well suited to symmetric, non-negative definite matrices. The PSEUP entry 
assumes this. Matrices formed by computer arithmetic will not always be symmetric. Hence, the 
routine always forces the symmetric matrix B to actually be symmetric, by taking the average of 
the element and its transpose. The nonsymmetric entry, unfortunately, approximately squares the 
ratio between largest and smallest eigenvalue. There is a nonsymmetric feature. The routine choses 
AA * , instead of the other way around, if A is a square matrix. This arbitrary choice agrees with the 
DECOM routine and the ASP routines. As a result, in the singular case, multiplying A by its pseudo- 
inverse from the left is more likely to give a diagonal matrix of 1 ’s and 0’s, than multiplying from 
the right side of A. 

Pivot Size 

DEP1 is used to compute a “smallest allowable pivot.” In no case is it reasonable or desirable 
to worry about exact equality in the use of such tolerances. Fortunately, work with ill-conditioned 
systems shows a series of pivots that decrease steadily in magnitude. Furthermore, the first “bad,” 
erroneous pivot is, at most, 10 to 1000 times smaller than its predecessors. Since ANDRA is choos- 
ing largest pivots first, the user has considerable latitude in actual choice. All positive elements can 
be accepted, if the matrix is known to be nonsingular, by choosing DEP1 very small. 

By choosing DEP1 very large — say, nearly 1.0 — the routine can be forced to reject pivots after 
the first. At present, there is no way of making it start iterating without having found at least one 
pivot. In other words, ANDRA always finds all the pivots it can before any side calculations are done. 
If this first rank is maximal, it never iterates. The first pivots are not in doubt, so these rules are more 
efficient. The routine always uses a tolerance for pivot acceptance; however, it uses a new tolerance 
50,000 times smaller than the last pivot found, for each call to find one pivot in iterative mode. The 
expensive test of matrix norms is avoided when no new pivot occurs. The PSEU routine has only a 
finite number of tries to find a new pivot before it quits. The exact number is the same as the maxi- 
mum rank. Since ANDRA has usually found several pivots initially, this is ample. 
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Iteration 


If DEP2 is larger than 1, it is raised to a power, used as a factor, and tends to make the 
routine stop with a smaller rank. DEP2 of 1 actually works for most iterations. 

Subroutine AN DR A 

The basic algorithm can be used as a separate routine by itself (see ANDRA documentation). 
The routine requires considerable setting and testing of parameters. It has an escape exit for too 
many iterations (calls to find only one pivot) without finding any. It returns a matrix, T, such 
that, if X if pseudoinverse of positive definite matrix A, then 

T l T = X 


Accuracy 

In double precision, the accuracy has been very good. Maximum accuracy can be obtained by 
using symmetric matrices and the PSEUP entry. The test program included in this manual as 

example 6 shows errors (determined by calculating AA&A—A# and A#AA#—A) on the order of 
1 0" 1 5 or less. 

The routine was also tested on the ill-conditioned 7X7 matrix in the ASP manual (NASA 
CR 475, p. 1 50). The exact inverse is given on page 151, and the error obtained from the ASP 
program using the equivalent of the PSEUP entry (p. 152) is on the order of 10' 1 . The error 
obtained using the VASP program and the PSEU entry was on the order of 1 0‘ s or less. 

Singular Case 

The routine forms a new inverse from the original symmetric matrix. Since there are several 
steps more between the inverse and the original input A, it is only natural that accuracy should fall 
off. In many cases, this inverse will give a diagonal matrix of O’s and l’s when used as a left inverse 
of A (or possibly as a right inverse). The work of reinverting B requires no extra matrices; it does 
destroy the usual values of C and D. No iteration can be done in the stage after B is found to be 
singular. It can be asked for in the starting stage. Error exits are made if the rank changes during 
reinversion. The smallest allowable pivot is redetermined. 

Error Exits ( Messages ) 

The error exits are reasonably self-explanatory. Unless otherwise noted, the errors cause a 
return from PSEU without completion of the calculations. Subsequent calculations in other portions 
of the program are suspect. 

Message Explanation 

Dimension error The total number of matrix elements was too 

large or too small. 


74 


Diagonal elements of matrix = 0 


Symmetric matrix B has no positive diagonal 
elements. Check input A. 


Rank has decreased 


Rank has increased 


Rank greater than matrix size 


Singular case. Reinverting, and the new rank is 
less than that of the earlier phase. 

Singular case. Reinverting, and the new rank is 
greater than in the earlier phase. Computation 
continues. 

RANK returned from ANDRA is greater than 
maximal rank. 


Timing 

The ANDRA routine by itself is very fast. The iteration mode is slower by a large factor than 
the regular mode of subroutine PSEU. 

The time estimates below (in hundredths of a second) are as used on the NASA Ames 360/50. 
High and low estimates are given, because real-time figures reflect an unknown percentage of time 
devoted to another CPU user. 


Case 

High 

Low 

PSEU, 2X2 matrix 

2 

i 

PSEUP, 4X4 matrix, reinvert 

14 

10 

PSEU, 7 X 7, no pivot rejection 

42 

30 

PSEU, 7X7, rank 6, reinversion 

103 

62 

PSEU, 7X7, iteration, no tests 

53 

30 

PSEU, 7X7, iteration, one test 

182 

118 

PSEU, 7X7, iteration, some tests of pivots 

253 

170 

PSEU, 7X7, iteration, many pivot tests 

501 

286 

PSEU, 7X7, iteration, nearly all tests 

607 

324 

PSEU, 4X2, reinversion 

3 

1 

PSEU, 4X2, reinversion 

2 

2 


METHOD 

Summary of Method 

PSEU has two entry points. The nonsymmetric entry forms A*A or AA t , whichever is 
smaller. At the end, A* is used again to form the pseudoinverse. Square A uses A*A. The result 
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is always forced symmetric afterward, even for symmetric entry. ANDRA is called to diagonalize 
this result in B. Most of the pivots are found and the steps made on the first call. If not iterating, 
this part is not repeated. If singular (rank of symmetric input not maximal), a transforming matrix 
is computed. A copy of the original symmetric B is transformed and reinverted by ANDRA. The 
result is retransformed by premultiplication and postmultiplication. If iterating, the pivot tolerance 
is decreased and ANDRA is called to find one pivot at a time. A side calculation is done to measure 
the quality of pseudoinverse formed at each step. The routine backs up one step and stops with 
rank one less if it makes a bad step. The result, if singular, is sent through the reinversion above. The 
use of PSEUP by DECOM avoids reinverting in the singular case, also it never uses a nonsquare input. 
There is a “find rank only” option. 

If PSEU is used without iteration, four I/O matrices are needed plus a dummy matrix. 

Iteration uses four additional dummy matrices. Iteration cannot be done during the reinversion. 
Besides those mentioned, entries BDNRM, MULT, and NORM are used for iteration. TTRM is also 
used except in rank only case. 

ANDRA (diagonalization algorithm). For a detailed description of the method, see the 
documentation of ANDRA itself. A mathematical description and examples are given in NASA 
CR-475. Subroutine PSEU calls ANDRA to do each pivoting step, after first forming a symmetric 
matrix B, which is indeed forced to be perfectly symmetric. 

The first call of ANDRA is an initialization call. An identity matrix T is formed. The rank 
counter is set to zero. On an initialization call, the routine proceeds to search the diagonal for 
pivots, as always. But after finding a pivot, it always goes back and looks for another pivot, 
regardless of the iteration option. The process of searching for pivots continues until the number of 
tries is one greater than the row size (no such test is made in the iteration case). If the routine fails 
to find a single pivot in the initialization call, it exits with an error message. Pivots are accepted if 
and only if they are not less than a threshold input at every call. Supposing that a pivot has been 
found in the diagonal, the next step is always the same. First the pivot is reduced to unity. That is, 
both the pivot row and column are divided by the square root of the pivot in B. Only the row of T 
is so reduced. The next step is to eliminate the pivot coefficient from all other rows not yet used as 
pivots. This part is the same as in other inversion methods. Both B and T are treated exactly alike 
here. Note that the actual algorithm checks the diagonal of a row to see if it is already marked as a 
pivot. If so, that entire row, and the row in T, are skipped. The pivot is then marked by an artifi- 
cial code. The routine always tests for this code and does not use this row again. The code is put in 
the actual pivot position. Thus the rows and columns are left in their starting places in the working 
matrix B. PSEU converts the result to a matrix of 1 ’s and 0’s that shows the independent and 
dependent variables. 

The code is tested for an integer. This is a considerable economy. The resulting T is never 
singular. If B were nonsingular and X the desired inverse of B, 

T*T = X 
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This part is done by subroutine entry TTRM, using coding shared with the iteration method. The 
final answer is put back in matrix B. (PSEU always uses the original A again rather than the origi- 
nal B, after this to give an answer for A. Thus, ANDRA is always supplied with a symmetric 
matrix B.) 

If B were singular at the start, a further reinversion would have to be done. See the next 
section. 

The Singular Case 

Suppose that the rank of B in the diagonalization by ANDRA does not turn out to be maximal, 
then PSEU must perform a number of matrix multiplications and call ANDRA and TTRM to reinvert. 
The accuracy is bound to suffer, even though the reinversion is done on an exact copy of the original 
B. A very short justification is given belov , followed by a close description of how the work is 
actually done. 

There exists a permutation matrix P, such that 

E = PTBT^ 1 

is a matrix of 0’s and 1 ’s (were it not for round-off error), with all the 1 ’s contiguous, starting in the 
first diagonal. If B had been so permuted before diagonalizing, then this different T resulting 
would be the one that gives an inverse that corresponds correctly to the old. But, since one is using 
a premultiplication and a postmultiplication, simple substitution of a permuted matrix does not 
work. (It would if matrix multiplication were commutative.) Thus, if it is necessary to transform 
the original starting B, reinvert, then transform back again. 

The permuted form of T (which does not actually occur) has a nonsingular corner submatrix, 
followed by the rest of diagonal set to l’s. These latter l’s correspond to the dependent equations 
of the original. 

The rule for constructing the transforming matrix U is given below. This matrix is made from 
T and put into the same storage T. The explicit construction of U is more efficient (in FORTRAN). 
From here on, the explanation concerns what is actually done, rather than the mathematical reasons. 

Let dj denote the ith diagonal element of B. (In case the reader has forgotten, this has been 
changed to a diagonal matrix of 0’s and l’s.) Given T, there are two cases: 

Case One: For Ujj not on the diagonal 

Use -tjj, if dj = 0; 

Use 0, if dj = 1 

Case Two: For Ujj on the diagonal 

Use the corresponding value of dj 

Next, using a copy of the original B, form 


C = U*BU 
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The result is actually put in the same storage that held B originally. The smallest allowable pivot 
for ANDRA is recalculated. This result, C, is sent to ANDRA to do the diagonalization again. The 
fact that C has rows and columns of 0’s that ANDRA has to skip makes the diagonalization ineffi- 
cient, but this cannot be helped. No iteration is done here. Let T 2 denote the result of this second 
ANDRA call. Then the new pseudoinverse is: 


X 2 =T 2 l T 2 


Transform this back to get a correct answer: 

X-UXaU 1 

The rest of the computation is as usual. Note that if the rank changes in the second ANDRA call, 
error exits are taken. 

Iteration 

The main method itself is purely algebraic. The iteration option is a way of estimating the 
amount of error in the generalized inverse and using this to stop with a smaller effective rank. Let 
B denote a matrix and X its pseudoinverse (after taking so many pivot steps in ANDRA). Then the 
two Moore-Penrose axioms read: 


BXB = B 
XBX= X 

If the iteration mode is selected, ANDRA first finds all the pivots it can. Then subroutine BDNRM 
is called twice. Each call returns the value 

norm(Q*P*Q - Q)/norm(Q) 

The values of P and Q are B and X in one call, X and B in the other. The resulting two small 
scalars (which would be zero if the axioms were perfectly satisfied) are added together. The result 
is taken as a factor times DEP2 raised to the current number of pivots. From successive iterations, 
one obtains a sequence of positive numbers, decreasing as one approaches the largest possible rank. 
As long as the new result is not larger, then a new pivot is searched for. If not, PSEU reverts to the 
previous values, before the current pivot was used. 

In practice a number of modifications are made. First, the pivot used last is returned as DEP3, 
even if rejected, so that the user can reconsider acceptance of it. Second, if maximum rank is 
achieved prior to iteration, no side calculations are done. Third, the smallest possible pivot allow- 
able is set to 0.00002 times the most recent pivot in order to reject many spurious pivots without 
doing the lengthy side calculations. This modification is based on actual observation of pivot 
behavior. The successive pivots of an ill-conditioned matrix usually decrease fairly rapidly. But 
there is usually a hugh jump in order of magnitude between the last good pivot and the first bad one. 
Parts of the side calculations are actually done in single-precision, to save time. Please note that a 
single iteration, besides the ANDRA call, makes ten subroutine calls, and one library routine call. 
Naturally, this is slow. 
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Matrix Storage Flow 

This section uses the same names as the Fortran IV routines. It tells what is put into each 
matrix of PSEU at various times. The call is CALL PSEU(A,B,C,EE,DEP,IP,D). The matrices A 
and B are the same size (possibly nonsquare). Matrix C is square with dimension equal to the 
smaller dimension of A. The other matrices are the same size as C. Matrix D is divided into five 
matrices. Let these be denoted as D, Dl, D2, D3, and D4. The last four are used only in iteration. 

Maximal Rank Case 

A symmetric matrix from A is placed in B (either directly, as in PSEUP, or indirectly, from 
matrix multiplication). A copy of B is put in D, unless the rank only, no iteration is used. ANDRA 
is called to diagonalize B and place the result in C. 

If the result is accepted, TTRM puts the generalized inverse of B into EE. Then the inverse 
of A is put into B. The A transpose may have to be used to get an answer for A. 

Singular Case 

The matrix U of the method is computed from C and put into C. (D holds original B.) 

EE = C 1 X D 
B = EE X C 

ANDRA is called to diagonalize B. Answer goes to EE. TTRM puts pseudoinverse of B into D. 

B = C X D 
EE = B X C l 

The pseudoinverse is now in EE, where the maximal rank case puts it. Routine now forms 
pseudoinverse of A in B. 

Iteration 

Before each call of ANDRA the current values of B and C are stored in Dl and D2, 
respectively. B and C are changed when a new pivot is used in ANDRA. BDNRM computes a 
number to decide if the pivot is to be rejected. EE, D3, and D4 are used as working storage in 
BDNRM. EE actually has a matrix put in it that would be zero if the Moore-Penrose axiom were 
perfectly satisfied. If the pivot is rejected, the old values from Dl and D2 are put back into B and 
C. The work of the singular case is done next if the call was not made from DECOM. 

Rank Only 

If iteration is used, a full complement of matrices must be used. In the ordinary case, 
matrix D may be omitted, and also matric E is not used. Naturally, no pseudoinverse is returned. 
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Figure 7.- Information Systems Co. flow chart - subroutine PSEU (A,B,C,EE,DEP,IP,D) 
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Figure 7.- Concluded 























29. BDNRM 


DESCRIPTION 

This subroutine computes the quantity 

norm(QPQ^ — Q)/norm(Q) 

where the values of P and Q are in the square arrays CT and EE or EE and CT, depending on the 
sign of NR. If P = Q#, the return value is zero. This routine can thus be used to test the quality of 
a pseudoinverse. 

USAGE 

CALL BDNRM(NR,CT,EE,D,KRV) 

Input Arguments 

Matrices: CT, EE with dimensions NR X NR 

Constants: NR, size of matrices and the sign controls multiplication procedure 

Output Arguments 

None: This is a function subroutine 

Dummy Arguments 

Matrix: D dummy array of size 5*NR 2 

Constant Array: KRV designates location of submatrices of D 

KRV1 = NR 2 
KRV 2 = 2*NR 2 
KRV3 = 3*NR 2 
KRV4 - 4*NR 2 


30. ANDRA 
SUMMARY 

ANDRA is a Fortran routine to diagonalize a positive definite symmetric matrix. The routine 
was originally designed to be used by subroutine PSEU. The routine has a parameter to command 
it to initialize on the first call. Two different modes can be used for pivoting steps. In the first 
mode, the routine does only one pivot to eliminate only one row at a time. In the second mode, as 
many pivots as possible are done in one call. Pivots are chosen in order of decreasing magnitude. 
They are rejected if smaller than a parameter threshold. The original matrix input is destroyed and 
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replaced with artificial values. However, symmetry is kept after each pivot. The answer matrix, T, 
is such that if X is the inverse of the input, 


X = T t T 

The routine has error exits for matrices of the wrong size, and for those that allow no pivot on the 
first try. 

USAGE 


CALL ANDRA(B ,T,DPR, JP) 

Name 

B 

T 

DPR 

DPR1 

DPR2 

JP 

JP1 

JP2 

JP3 

JP4 

JP5 


Description 

Input symmetric matrix. Destroyed. 

Answer. T*T = inverse of B. 

Parameter array of size 2. 

DPR1 is the tolerance for trial pivots. Any less 
than this are rejected as zero. 

DPR2 is the last pivot actually used. Unchanged 
if no new pivot found. 

Integer parameter array of size 5. 

Zero if all pivoting to be done on one call; 
nonzero if only one pivot per call. 

Zero if initialization call. Subroutine sets to one 
when a pivot is found. 

Holds the effective rank = number of pivots 
found. 

The integer giving the row and column size. May 
range from one to a nominal figure. 

The integer row where the last pivot was found. 
The rows are left in the same positions as in the 
input matrix. 
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Figure 8.- Concluded. 



















METHOD 


Mathematical 

The method is described in the ASP manual, pages 137 to 139. 

The square matrix T is initialized to be the identity. 

Step 1 

The diagonal of B is scanned to find the largest pivot. Pivots are only taken from the 
diagonal. If no pivot is found, skip to step 3. 

The square root of the pivot is taken. The pivot row and pivot column are divided by the 
square root. Thus, the pivot, at the intersection of the row and the column, is reduced to unity. 
The corresponding row of T is also divided by the square root. 

Step 2 

The new, reduced pivot row is used to eliminate the elements of the pivot column. Let K be 
the pivot row and column. The pivot row is multiplied times the element in the j,k position. The 
resulting row-vector is subtracted from the jth row. This process is repeated for each row j that 
has not yet been a pivot row. Exactly the same operations are carried out on the corresponding 
rows and columns of T, except that, of course, the multiplier for a pivot row comes from B. Then 
the pivot row of B, except for the pivot, is set to zero. The pivot row and its corresponding row in 
T are never used again. 

Step 3 

If the rank is maximal, exit. If no pivot has been found, a test is made to see if this should be 
an error exit, or normal exit. Otherwise, repeat step 1 . 

Computational 

In practice, a number of modifications are made. The actual calculations are rewritten to 
optimize speed and storage. The reciprocal of the square root is used, instead of a division. For 
single precision, straight division would probably be best. In step 3, an artificial code is put into the 
pivot position. This code is chosen as one that cannot be the result of floating-point arithmetic. 
Such a technique works in a great many different Fortrans. If a row is found to be marked by a 
pivot code, it is skipped in steps 1 , 2, and 3 above. 

The pivot position is forced to be exactly 1 before step 2 is applied. The pivot-code is actually 
tested for as an integer. The pivot size is tested for in single precision. These modifications are for 
speed. A count is kept of the number of pivot searches. If this count is one greater than the num- 
ber of rows, the routine always stops searching for pivots. The result, if B has maximum rank, is 
a matrix T such that T*T = inverse of B. The input B consists of 0’s everywhere except the 
diagonal, which holds pivot codes. 
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Error Messages From ANDRA 
Message 

Dimension error 


Finds no pivots 


Explanation 

The total number of matrix elements was too 
large or too small. The parameter JP(4) cannot 
be less than one nor more than MAXRC. 

ANDRA could find not a single pivot in its 
very first search of diagonal. 


31. DECOM 
SUMMARY 

Fortran IV subroutine DECOM generates four double-precision output matrices from the 
symmetric, non-negative definite input matrix B. One output is a matrix S such that if E is a 
unity matrix of rank r, then 


B = SEE^ 

This matrix is obtained as the inverse of a matrix T, by calling subroutine INV; T comes from 
subroutine PSEU. It is defined by TBT* = E, a diagonal matrix with elements 0 or 1 . E is also 
returned, along with a permutation matrix P such that 

PEP 1 = I r 

a diagonal matrix with all 1 ’s moved to the uppermost left corner. Given these matrices, and the 
ability to find a pesudoinverse of A, a decomposition of any matrix is possible. PSEU and DECOM 
are called and the matrices then multiplied as described in the method to give a Kronecker decom- 
position. The routine calls PSEUP and INV to do most of the calculation. Besides returning the 
matrices P and E, it does nothing that could not be done by successive calls of other matrix 
routines. It has parameters and error exits similar to that of PSEUP. 

USAGE 

CALL DECOM(A,B,C,E,J,DCM,KP,D) 

Arguments Description 

A The symmetric non-negative definite input. 

B The output matrix E, diagonal of 0 and 1 , with 

l’s in the independent columns. B, C, E, J, D, 
and D 1 are all of same size as A. 

C The output T, such that TAT* = diagonal of 

0’s and 1 ’s. 
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E 

J 

DCM 

DCM1 

DCM 2 

DCM 3 
KP 
KP1 
KP2 


KP3 

KP4 

D 


Holds the inverse of A (B does not hold the 
inverse of A). (Not E of ASP.) 

A square integer matrix for housekeeping in 
INV and DECOM. 

Parameters, just' as in subroutine PSEUP. 

Multiplied times the largest magnitude of diagonal 
of A, to give a lower limit for an acceptable pivot 
in PSEUP. Set at 2( 1 0) -6 if zero is input. 

Used only if the user elects to iterate in PSEUP. 
Set at 1 .DO (no effect) if zero is input. 

Note: DECGEN uses the default options for 
DCM1 and DCM2. 

The last pivot accepted by subroutine PSEUP, 
during diagonalization of input matrix A. 

Integer control parameters, just as for subroutine 
PSEUP. 

Zero, do not iterate in PSEUP. One, iterate in 
PSEUP. 

Zero, do all calculations. Nonzero, do rank only. 
Changed to reflect the rank on output. Should 
be set to zero or minus one before each call. 

Note: DECGEN uses KP1 and KP2 = 0. 

The row size of the matrix input. 

The column size. 

Note: This parameter is forced negative as a 
signal if T cannot be inverted by INV. 

D has five parts, as does the “dummy” array in 
PSEUP. Let these be denoted D, D1,D2, D3, 
and D4. These five equal arrays must be included 
in the size of parameter D if iteration by 
PSEUP is selected. If no iteration is used, D2,D3, 
and D4 may be omitted. D holds the inverse of 
output C. D1 holds the permutation matrix P. 
Note: If rank only is computed, D1 is computed, 
but D is not. A, B, C, and D1 are thus the only 
matrices with useful values returned. 



METHOD 


The results from DECOM are an effective rank r; matrices B and D, which are used in 
further calculations to get a Kronecker decomposition, or to see which variables are dependent; and 
the permutation matrix P inDl. This section describes the sequence to obtain the Kronecker 
decomposition in two different cases. The goal is two matrices G and H. DECOM does not produce 
these matrices; they are produced either by DECGEN or by the user according to the following steps. 

Let R be the matrix to decompose. Matrices G and H are desired such that 


R = HG 


H is to have only r nonzero columns; G is to have only r nonzero rows. Small r is the rank of 
R. 

Case 1 

Matrix R is symmetric, non-negative definite. Input R as the square input A to DECOM. 
Then H and G are produced afterward from the matrices in the call statement as follows: 

Parameter B is a diagonal matrix with r l’s; H and G are computed by : 

H = D X B 
G = (D X B)* 

A = original = DBB*D* 

Case 2 

R is nonsymmetric, possibly not even positive definite. Form RR* (subcase a) or else form 
R*R (subcase b). The subcases are chosen to give the smaller dimensions. If R is square, use RR l 
to agree with both PSEU and DECOM. Let this symmetric result be the input A to DECOM as 
in case 1. Obtain D and B as before and save them. In subcase a,X=R t XE, but in subcase b, 

X = E X Rt. Then for subcase a, take 


H = DB 
G = (XDB) 1 


Similarly, in subcase b, take 


H = X f DB 
G = (DB) 1 

Note: The H and G matrices produced have the same dimensions as the smaller dimension of R. 
If the rank of R is not maximal, there will be zero rows or columns in H and G. If the matrix D1 
is used instead of B in the above calculations, the zero rows or columns will be at the right or bot- 
tom, and the dimensions may be easily reduced. This latter is the procedure used in subroutine 
DECGEN. 
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Computation 

In practice, the subroutine is very short; it calls on PSEUP and INV to do the computations. 
No flowchart is needed, since there are no loops of any consequence. 

Step 1 

The matrix size is tested for reasonableness, with an error exit if it is not. KP(1) is set to 
special negative values to suppress reinversion by PSEUP, and to change somewhat the matrix 
outputs. This change is not discussed in PSEUP. 

Step 2 

Entry point PSEUP is used to diagonalize the input. C holds a matrix T such that 
TAT* = B, a matrix of 0’s and l’s. If the rank only option is input, the routine skips to step 4. 

Step 3 

Subroutine INV puts an inverse of T into D. The flag PIV is tested. If zero, INV 
failed; the routine prints an error message. INV uses matrix J. 

Step 4 

The matrix E, which is matrix of 0’s and l’s, is scanned along its diagonal. A matrix P of 
0’s and l’s is constructed such that 


PEP 1 = I r 

I r has all Fs moved to extreme upper left corner. A record of successive diagonal positions that 
are 0 is kept. As each 1 is found in the diagonal in position k, the record is checked to see if 
there is an earlier 0 (or 1) that needs to have a 1 permuted into its place j by permutation p. 

If so, a 1 is put into position j, k of P. Premultiplication by P will move position k, k to j, k. 
Postmultiplication by P* will move j, k to position j, j. Position k, k is also marked as a 
hole that could be filled by a 1 lower on the diagonal, since it vacates its old position. The 
record in the first column of J has the structure of a queue. Matrix P is in D 1 , the second 
matrix of dummy array D. 

Step 5 

Return. 

NOMENCLATURE 

The nomenclature used in DECOM is compatible with that used in PSEU, but differs 
from that used in the ASP manual description of the decomposition routine, p. 154. Also, 
since DECGEN requires dummy storage, the nomenclature there is different again. The 
following table lists the correlations: 
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I 


DECOM 

DECGEN 

ASP 

A 

DUM1 

aa t 

B 

DUM(N7) 

E 

C 

DUM(N4) 


E 

DUM(N5) 


D 

DUM(N2) 

S 

D1 

DUM(N3) 

P 

J 

DUM(N6) 




APPENDIX B 


LISTINGS OF ALL VASP SUBROUTINES 



SIJBR PUT IMF READ (I, A, NA, B , NB, C, NCX, 0, NO , E, ME) 


DIMENSION All), B(l), C(l), 0(1), E ( 1 ) 

D I M F M S I ON N A ( 2 ) , N B ( 2 ) , M C X ( 2 ) , Nl)( 2 ) , NE ( 2 ) , MZ ( 2 ) 
DOUBLE PRECISION A, B, C, 0, E 

READ ( 5, 100 ) 

LAR, 

NZ ( 1 ) , 

NZ ( 2) 

CALL READKA 

, MA,MZ, LAB) 



I F ( I .EO. 1) 

GO TO 999 



READ! 5,100) 

LAR, 

NZ ( 1 ) , 

NZ (2 ) 

CALL RE ADI ( B 

, NB,MZ, LAB) 



I F ( I .FO. 2) 

GO TO 999 



READ ( 5 , 100 ) 

LAB , 

NZ ( 1 ), 

NZ ( 2) 

CALL RFADKC 

, MCX ,NZ , LAB) 



I F ( I .EO. 3) 

GO TO 999 



READ! 5,100) 

LAB, 

MZ ( 1) , 

NZ ( 2 ) 


CALL REA01 (D, MD,NZ, LAB) 

TF ( T .FO. U GO TO 999 

RFAO ( 5, 100 ) LAB , MZ( 1 ) , NZ ( 2) 

CALL RFADKF, NF t NZ, LAB) 

TOO FORMAT! A4.4X.2T4) 

99 9 RETURN 
FNn 


SUB R OUT I NF RDTITL 

COMMON /LINES /NLP, LIN, TITLE (23) 
RE AO (5,100) (TITLE ( T ), 1 = 1, IK) 

100 FORMAT (18A4) 

CALL LNCNT(IOO) 

RETURN 

END 


VO 

u> 



SUBROUTINE PRNT ( AR.NAR ,NAM, IP) 

C SUBR PRNT PRINTS ROUBLE PRECISION MATRIX 
COMMON /FORM/NEPR,FMTl(6), FMT2 ( 6) 

COMMON/ LINES/ NLP. LIN* TITLE! 23 ) 

COMMON /MAX/MAXRC 

C- NOTE NLP MO. LINES/PAGE VARIES WITH THE INSTALLATION. 

DATA KZ , KW , KB /1HQ,1H1,1H / 

REALMS AR 0 

DIMENSION A R ( 1 ) , NAR ( 2 ) 

NAM F = MAM 1 

C-IF IP =1, HEADLINE SAME PAGE, IF IP =2, HEADLINE, NEW PAGE 

C I P = 3 , NO HFADLINE, SAME PAGE, IP = 4, NO HEADLI NE» "NEW PAGE 

II = IP 

NR=NAR( 1) 

MC=N AR ( 2 ) 

NLST = NR * NC 

IFINLST .GT. MAXRC .OR. NLST .LT. 1 . OR . MR . L T . 1 ) GO TO 16 
I F ( NAME .EO. 0) NAME = KB 
C- SKIP HEADLINE IF REQUESTED. 

GO TO ( 11,1TTTT3ZTT'2) ; TT 

10 CALL LNCNTI 100) 

11 CALL LNCNT ( 2 ) 

3 WRITE (6, 177 ) KZ , NAME , NR , NC 

177 FORMAT! A'l ,5X,A4,8H MATR I X , 5X , I 3, 5H ROW S , 5X , I 3 , 8H COLUMNS) 

GO TO 13 

12 CALL LNCNT ( 100 ) 

GO TO 13 

132 CALL L NCNT ( 2) 

WRITE (6,891) 

891 FORMAT (1H0) 

C- BELOW COMPUTE NR QF LINES/ ROW 

13 J=(NC-1) /NEPR+1 

NLPW = j 
J ST = 1 

C- COMPUTE LAST ROW POSITION -1 

NLST = NLST -MR 

MN=NC 

IF (NC.GT.MEPR ) MN = NE PR 
KLST = NR*( MN- 1 ) 



91 

CONTINUE 
DO 912 J = JST , 
CALL LNCNT(NLPW) 

NR 




KLST = KLST +1 
WRITE ( ft , F MT 1 ) 

( AR ( M ) , N = J, 

KLST, 

|M U ) 


IF (NC.LF.NFPR) 

GO TO P 1 2 




ML ST = NLST +1 
KMR=KLST+MR 





WRITE (6.FMT2 ) ( AR (N ) , N=«NR .NLST , 

NR ) 


912 

CONTINUE 

RETURN 




1ft 

CALL LNCNT(l) 





WRITE (ft, 91ft) 

MAM t MAR 



91ft 

FORMAT (' FRRPR 
CALL ASPERR 

I N P R NT MATRIX 

' A 4, ' 

HAS NA= ' , 21ft) 

RETURN 

FND 


SUBROUTINE LMCNT ( M ) 

COMM flM/L INFS/ ML P ,L IM,TT TLF ( 23) 

L I M = L I N+N 

IF ( L I M « L F . M L P ) GO Tn 2 0 

WRITE (ft, min) (TITLF(I),T = 1,23) 
1010 FORMAT (1H1,23A4/) 

L I N=2 + N 

IF (M.GT.MLP) L T im = 2 
20 RETURN 
FND 


i 


II BR GUT I N E ADD ( A ,N A ,B , NB , C , NC ) 

ImEMSIOM a { rTTBl n 7tr I J ,NA(2), NB( 7T7WT7 ) 


COMMON /MAX/MAXRC 

DOUBLE PRECISION A,B,C 

IF( (NA( 1) ,NE.NH(1 ) ) . OR. (MAC 2) .ME. NR ( 2 ) ) ) GO TO 999 
NC( 1 )=NA( 1 ) 

NC ( 2 ) =NA ( 2 ) 

L=NA ( 1 )*NA ( 2 ) 

IF ( NA( 1) .LT . l.OR.L.LT .l.OR.L.GT . MAX RC ) GO TO 999 
DO 300 1=1, L 

30"0 C( n=A(T)+B(I ) 

GO TO 1000 
999 CALL LMCNT (1) 

W R I T E ( 6 * 50 ) NA,NR 

50 FORMAT (' DIMENSION ERROR IN Al >o N A = • 2 I 6 , 5 X, ' NR= ' 2 I 6 ) 

CALL ASPERR 


1000 RETURN 
END 



SUBROUTINE SUBT(A f NA,B,NR,C, NC) 

DIMENSION A(I) t B(l),C(l),MA(2)tMB(2) ,NC(?) 

DOUBLE PRECISION A,B,C 
COMMON /MAX/MAXRC 

I F ( ( N A ( 1 ) , NE « NB ( 1 ) ) . OR » ( M A ( 2 ) »NE » MR ( 2 ) ) ) 00 TO 999 

NC( 1 )=MA( 1 ) 

NC ( 2 ) =NA ( 2 ) 

L =NA ( 1 ) *N A { 2 ) 

TF — rMA m . LT . 1 . OR . L . L I . 1THR : OrTTTOTTTJ CM II) 999 

DO 300 I =1 » L 

300 C ( I ) = A ( I ) -R ( I ) 

GO TO 1000 
999 CALL LNCNT ( 1 ) 

WRITE(6«50 ) NA,NB 

50 FORMAT (• DIMENSION ERROR IN SUBT NA = • 2 16 r 5 X , ' NB= ' 2 1 6 ) 
CALL ASPERR 

1000 RFTIIRN 

END 



'JD 

00 


SUBROUTINE MUL T ( A ,NA , B ,NR , C , MC ) 

DIMENSION A!1),B(1),C(1),NAI2),NB<2) ,NC(?) 

DOUBLE PRECISION A,B,C 
COMMON /MAX/MAXRC 
NCI I ) = NA ( 1 ) 

NCI 2) = MRI2) 

IFINAI 2) .NE.NB(l) ) GO TO 999 
NAR = MA(l) 

NAC = NAI2) 

NBC = NBI2) 

NAA=NAR*NAC 

nbb=nar*nbc 

IF ( NAR . LT. 1 . OR. NAA.LT.l. OR. NAA.GT.MAXRC.OR.NBB.LT. l.OR. 

1 NBB.GT.MAXRC) GO TO 999 
IR = 0 

I K=-NAC 

DO 300 K= 1 » NBC 
IK = IK + MAC 

DJ 30 Q J= 1 tMAR 

I R=I R+l 

I B= I K 
J I = J -NAR 
C 1 1 R ) =0 . 

DO 300 1=1, NAC 

JI = JI + NAR 

I B= I B+l 

C 1 1 R ) =C I IR)+A( JI )*6( IB ) 

300 CONTINUE 

GO TO 1000 
999 CALL LNCNT (1) 

WRITE! 6,500) I NA I I ) , I = 1 , 2 ) , I MB I I ) , I = 1 , 2 ) 

500 FORMAT (' DIMENSION ERROR IN MULT NA = "' 2 I bX,' « NB= 1 2 I 6 ) 
CALL ASPERR 

1000 RETURN 

END 



JUIJn* LA./ I1IMC OL>Ml_U \ M y PM My My ' 'J u y M / 

DIMENSION A(1),B(1 ),NA(2),NB(2) 
COMMON /MAX/MAXRC 
DOUBLE PRECISION A, R, S 

NB( 1 ) = N A { I ) 

NB ( 2 ) =NA ( 2 ) 

L = NA( I )*NA(2) 


IF ( N A ( 1) .LT. 1.0R.L.LT. 1.0R.L.GT.MAXRC ) 
DO 300 1=1, L 

on TO 999 



300 

B ( I ) =A( I) *S 




1000 

RETURN 




999 

CALL LNCNT(l) 
WRITE (6,50) NA 




50 

FORMAT (' DIMENSION ERROR IN SCALE IM A = 

CALL ASPERR 

RETURN 

>216) 




END 





VO 

M3 



o 

o 


SUBROUTINE TRAMP ( A ,NA, B , MR ) 

DIMENSION A (1 ) , B (1) ,NA( 2 ) ,NB (2) 

DOUBLE PRECISION A,B 

COMMON /MAX/MAXRC 

NB( 1 )=NA( 2) 

NB( 2 ) =N A ( 1) 

NR=NAI 1) 

NC =N A ( 2 ) 

L=NR*NC 

IF (NR ,LT.l»OR»L»LT. l.QR.L.GT»MAXRC) CO TO 999 
IR = 0 

DO 300 1=1, MR 

I J= I-NR 

DO 300 J=1,MC 
I J=I J+NR 

I R = I R + l 

300 BUR )=A( I J) 

RETURN 

999 CALL LNCNT(l) 

WRITE (6,50) NA 

50 FORMAT (' DIMENSION ERROR IN TRAMP NA='2I6) 

CALL ASPERR 

RETURN 

END 



SUBROUTINE INV (A , NA , DET,L ) 

DIMENSION A (1 ), L( 1) ,NA(2) 

DOURLE PRECISION A, RET, RIGA , HOLD 
COMMON /MAX/MAXRC 

IF ( NA( 1 ) .NE.NA (2 ) ) GO TO 600 

C SEARCH FOR LARGEST ELEMENT 

DET= 1. 

N=NA ( 1 ) 

NSQ=N*N 

IF (N.LT.l.OR.NSO.GT .MAXRC) GU TO 600 

NK = - N 

DO 80 K= 1, N 
NK = NK + N 

L ( K ) = K 

NPK=N+K 
L ( NPK ) =K 

KK = NK + K 

RIGA = A ( KK ) 

DO 20 J= K, N 

I Z = N» ( J - 1 ) 

DO 20 1= K, N 
IJ = IZ + I 

10 IF( DABS( RIGA) - DABS ( A ( I J ) ) ) 15, 20, 20 

15 RIGA = A ( I J ) 

LIK) = I 

NPK=N+K 

L ( N PK ) = J 
20 CONTINUE 

X TMTFRCHAMGF ROWS 

J = L ( K ) 

I F ( J - K) 35, 35, 25 

2 5 KI = K - N 

DO 30 I = 1, N 
KI = KI + N 

HOLD = **A ( K I ) 

JI = KI - K + J 
A ( K I ) =* A ( J I ) 

30 A ( J I ) = HOLD 



102 


C INTERCHANGE COLUMNS 

35 NPK=M+K 

I -L (NPK) 

I F (I - K) 45, 45, 38 
38 JP = N*( I - 1 ) 

DO 40 J = 1, N 

JK = NK + J 
JI = JP + J 

HOLD = -A ( JK ) 

A ( JK ) = A ( J I ) 

40 A ( J T ) = HOLD 

C DIVIDE COLUMN RY MINUS PIVOKVALUE OF PIVOT EL FRENIS IS CONTAINED 

C IN RIGA) 

45 I F ( RIGA) 48, 46, 48 

46 PET = 0. 

CALL LNCNT (I) 

KKK=KK-1 

WRITE (6, 1046) KKK 

1046 FORMAT (' INV ERROR DETERMINANT OF A=0 RANK OF A=',J4) 

CALL ASPERR 

RETURN 

48 DO 55 1= 1, N 

IF (I - K) 50 , 55 , 50 

50 IK = NK + I 

AUK) =-A ( IK ) /( RIGA ) 

55 CONTINUE 

C REDUCE MATRIX 

DO 65 1= 1, N 
IK = NK + I 

HOLD E _A i_ LI<J 

IJ = I - N 

DO 65 J= 1, N 

TJ = TJ + N 

I F ( I - K ) 60, 65, 60 

60 I F ( J— K) 62, 65, 62 

62 KJ = IJ -I + K 

A ( I J ) = HOLD* A( KJ ) + A( I J) 

65 CONTINUE 

C DIVIDE ROW BY PIVOT 



KJ = K - M 
no 75 J= 1, M 

KJ = KJ + N 

I F ( J - K) 70, 75, 70 
70 A ( K J ) = A(KJ )/B IGA 

75 CONTINUE 

C PRODUCT OF PIVOTS 

DET =DET*B IGA 

C REPLACE PIVOT BY RECIPROCAL 

A CKK ) = 1 ./RIGA 
RO CONTINUE 

C FINAL ROW AND COLUMN INTERCHANGE 

K = N 

100 K = K - 1 

I F ( K ) 150, 150, 105 

105 I = L( K) 

IF (I - K ) 120, 120, 108 

108 JO = N*( K - 1 ) 

JR = N-MI- 1) 

DO 110 J= 1, N 
JK = JO + J 

HOLD = AUK) 

JI = JR + J 
AUK ) = - A ( JI ) 

110 A C JI ) = HOLD 
120 NPK=N+K 

J=L ( NPK ) 

I F ( J - K) 100 , 100, 12 5 
125 KI = K - M 

DO 130 1= 1, N 
KI = KI + M 
HOLD = AC K I ) 

JI = KI - K + J 

AC KI ) = - A ( JI ) 

130 A ( J T ) = HOLD 

GO TO 100 




104 


150 RETURN 

600 CALL LNCMT (1) 

WRTTF (6. 1600) MA 

1600 FORMAT (' INV REQUIRES SQUARE MATRIX NA=',?I4) 
CALL ASPERR 




105 


SU B R OUT I N F N ORM ( A , N A , AW OR M ) 
'fSTMFWS I ON MA I 2 ) ,A( 1) 


DOUBLE PRECISION A , ANORM, SUM , R OWi'iAX , COLFAX 
COMMON /MAX/MAXRC 


MA R = MA ( 1 ) 

NAC = MA(2) 

L = N AR* NAC 

IE (NAR .LT.l.OR.L .LT. 1, OR.L. GT.MAXRC) 0,0 TO 999 
COL MAX = 0. 

ROWMAX = n. 

FT ="() 

DO 300 I = 1 , NAC 
SUM = 0. 

DO 301 J = 1,NAR 
K = K + 1 

301 SUM = SUM + D AB S ( A ( K ) ) 

IF (COLMAX.LT .SUM ) COLMAX = SUM 

300 CONTINUE 

DO 302 I = 1,NAR 

SUM = 0. 

K = I - NAR 

DO 303 J = 1 , N AC 

K = K + NAR 

303 SUM = SUM + DABS ( A ( K) ) 

IF (ROWMAX. LT. SUM ) ROWMAX=SUM 

302 CONTINUE 

ANORM = DMIM1 (COLMAX, ROWMAX) 

RETURN 

999 CALL LMCNT (1) 

WRITE (6,50) MA 

50 FORMAT (' DIMENSION ERROR IN NORM IMA=»2I6) 

CALL ASPERR 

RETURN 

END 








o 

Os 


SUBROUTINE UNITY(A,NA) 

DIMENSION A (1),NA(2) 

DOUBLE PRECISION A 

IF(NA(1).NE.NA(2)) GO TO 999 

CALL SCALF ( A,NA ,A ,NA, 0. DO ) 

J = - N A ( 1 ) 

MAX = NA( 1 ) 

DO 300 1=1, MAX 
J = N AX +J+1 

300 A ( J ) = 1 . 

GO TO 1000 
999 CALL LNCNT ( I ) 

WRITE) 6 » 50) ( MA ( I ) ,1 = 1,2) 

50 FORMAT (• DIMENSION ERROR IN UNITY NA='2I6) 
CALL ASPERR 

1000 RETURN 


END 



SUBROUTINE TRCE ( A , N A , T R ) 

DOUBLE PRECISION A ( 1 ),TR 
0 I MENS I ON NA(2) 

COMMON /MAX/MAXRC 

IF (NA(1).NF.NA(2)) GO TO 600 
TR=O.DO 

N = NA ( 1 ) 

NN = N*N 

IF (M.LT.l.OR.NM.GT.MAXRC) GO TO 600 

DO 10 1=1, N 

M = I+N*( I- 1 ) 

10 TR=TR+A(M) 

RETURN 

600 CALL LMCNTd ) 

WRITE (6,1600) NA 

1600 FORMAT (» TRACE REQUIRES SQUARE MATRIX NA=',?I6) 
CALL ASPERR 
RETURN 
END 


o 

-j 





o 

oo 


SUBROUTTNF FQUA TE ( A . N A , B , MR ) 

DIMENSION A (1),B( 1 ) , NA ( 2 ) , NB ( 2 ) 

DOUBLE PRECISION A, B 

COMMON /MAX/MAXRC 

NB ( 1 ) = NA ( 1 ) 

NB( 2 ) =NA(2) 

L=N A I 1 ) *NA ( 2 ) 

IF (NA(l).LT.l.OR.L.LT.l.OR.L.GT. MAX R C ) 00 TO 999 

DO 300 1=1, L 

300 B ( I ) = A ( I ) 

1000 RETURN 
999 CALL LNCNT ( 1 ) 

WRITE (6,50) NA 

50 FORMAT (» DIMENSION ERROR IN EQUATE NA=«2I6) 


CALL ASPERR 

RETURN 

END 



SUBROUTINE J UXTC ( A , N A , B , MB , C , NC ) 

DIMENSION A(1),B(1 ),C( 1 ) , NA ( 2 ) ,NB ( 2 ) , NC ( 2 ) 
DOUBLE PRECISION A,B,C 

COMMON /MAX/MAXRC 

IF ( NA ( 1 ) . ME .NB ( 1 ) ) GO TO 600 
NC ( 1 ) = NA ( 1 ) 


NCI 2 ) =N A ( 2 ) +NB ( 2 ) 
L=NA( 1) *NA(2 ) 


NNC=NC ( 1 ) *NC ( 2 ) 


IF ( NA ( 1) »LT » l.OR.L.LT. l.OR.L.GT .MAXRC ) GO TO 600 
IF (NC(2).LT.l. OR. NNC.GT. MAXRC) GO TO 600 
MS=NA( 1 ) *NA ( 2 ) 


ML 


10 C ( I ) =A( I) 


ft 


1=1. MS 


MBS=NA( 1)*NB(2) 

DO 20 I=1 t MBS 

J=MS+I 

20 C ( J ) =B ( I) 

RETURN 

600 CALL LNCNT(l) 

WRITE (6,1600) NA, NB 

1600 FORMAT (* DIMENSION ERROR IN JIJXTC, NA= 1 , 2 16, 5X, 1 NB= 1 , 2 16) 
CALL ASPERR 
RETURN 


£Na 


o 

VO 



SUB ROUT IMF JUXTR! A ,NA, B.NB,C»NC) 

DIMENSION A(1),B(1),C(1),NA(2),MB(2),NC(2) 

DOURLE PRECISION A,B,C 

COMMON /MAX/ MAX RC 

IF(NA( 2) .NF.NRI2) )GO TO 600 
NC(2)=NA<2) 

NCI 1 ) = NA ( 1 ) +NB ( 1 ) 

L = N A ( 1 )*NA (2) 

MMC = NC( 1) *NC(2 ) 

IF (NA( 1 ) .LT.l.OR.L.LT. 1 . OR . L. GT. MAX RC ) GO TO 600 

IF (NC(2).LT.1.0R.NNC.GT.MAXRC) GO TO 600 
MCA=MA(2 ) 

M R A = N A ( 1) 

MR B= MB ( 1 ) 

MRC=MC( 1) 

00 10 I=l t MCA 

00 10 J= 1 » MR A 
K=J+MRA* ( 1-1 ) 

L = J + MRC*( 1-1 ) 

10 C ( L) = A ( K ) 

00 20 1=1 , MCA 
00 20 J=1,MRR 
' K=J+MRR*( 1-1 ) 

L=MR A + J + MRC 5 ' 5 ( 1-1 ) 

20 C(L)=B(K) 

RETURN 

600 CALL LNCNT(l) 

MRITE(6,1600) MA,NB 

1600 FORMAT!' DIMENSION ERROR IN JUXTR, MA= ' , 2 I 6 , 5X , ' NR= ' , 2 I 6 ) 
CALL ASPERR 

RETURN 

END 



SUBR OUT INF FAT ( A , hi A , TT , B, MB, C, NC, DUMMY, KOUM ) 

DIMENSION A ( 1 ),B( 1 ) , DUMMY) 1 ) , NA ( 2 ) , NB ( 2 ) , ND ( 2 ) ,C(1 ) ,NC( 2) 
DOUBLE PRECISION A, T, TT , ANAA, TMAX, B , DUMMY, C , S, SC 

COMMON /MAX/MAXRC 

NR = MA( 1) 

NCC = NA(2 ) 

NC ( 1 )=NR 

NC ( 2 ) =NCC 

NB ( 1 )=MR 

NB( 2 ) =NCC 

LD=NR*NCC 

IF (MR.ME.MCC.OR.NR.lt. 1 .UR.LD.GT .MAXRC ) GO TO 998 

NDD = ?*NA( 1 )*NA (1 ) 

IF(KDUM .LT.NDD) GO TO 998 
NDD= NA( 1 )*NA ( 1) + 1- 
T =TT 

CALL NORM ( A ,NA , ANAA ) 

TMAX= 10.01/ANAA 

K = 0 

101 IF ( TMAX-T ) 103,104,104 

103 K=K+ 1 

T = TT / 2**K 

IF (K- 1000 ) 10 1,102,102 

104 SC=T 

CALL SCALE(A,NA,A,NA,T) 

DO 401 1= 1, 2 
401 NB ( I ) = NA( I ) 

CALL U M I T Y ( B , NB ) 

CALL SCALE (B,NB,DUMMY( 1 ),NB,T) 

S = T/2. 

CALL SCALE ( A,NA, DUMMY (N DO ) ,NA,S) 

N = 35 
1 1 = 2 

CALL ADD ( DUMMY ( 1 ), N A, DU MMY ( NOD ) ,NA, DUMMY (MDD) ,NA) 

CALL ADD ( A , NA , B , NB , DUMMY ( 1 ) , NA ) 

CALL EOUATF ( A , N A , C » NC ) 

106 CALL MULT ( A , NA , C , MC , B , MB ) 




S = 1 . DO/ I I 

CALL SCALE(R,NR,C,NC,S ) 

S = T / ( I 1 + 1) 

CALL SCALF(C,NC,R ,NB ,S ) 

CALL ADD IB, NB, DUMMY ( NDn ), MB , DUMMY ( NDD ) ,MR) 

CALL AUDI C,NC, DUMMY! 1 ) , NC ,DUMM Y ,N C ) 

N=N— 1 

IF (M) 107,107,105 

105 11 = 11 + 1 

GD TO 106 

107 CALL EQUATE ( DUMMY ( 1 ) , MB, B, NR) 

IF (K) 109,108,212 

109 CALL LNCNT (1) 

WRITE ( 6,110) 

110 FORMAT (' ERROR IN EAT K 15 NEGATIVE') 

112 IF ( K — 1 ) 213,212,212 

213 K=1 

212 DO 111 J= 1 , K 

T = 2*T 

CALL EQUATE ( R , NB , DUMMY ( 1 ) , MB) 

CALL MULTI DUMMY ( 1 ) , NA, DUMMY ( MOD) , NA,C,NC) 

Call add l diimMy (nDdT, NC,T, NC /"DUMMY (Nl)nj , W. ) 

111 CALL MULTI DUMMY! 1) , MR, DUMMY I 1 ) , NB , B, NB ) 

TT = T 

10B CONTINUE 

CALL EQUATE ( DUMMY ( NDD ) ,NC,C,NC) 

s = npo/sc 

CALL SCAL F I A, NA , A , N A , S ) 

RETURN 

1 02 CALL LNCNT (1) 

WRITE (6,119) 

119 FORMAT (' ERROR IN EAT K=1000') 

CALL ASPERR 

RETURN 

998 CALL LNCNT (1) 

WRTTFI6, 50) NA . KDUM » N DP . 

50 FORMAT! 'DIM E NS I ON ER R OR IN EAT NA='2I6, ' K[)1IM= ' I 5 , 5X , 

1 'KDUMI MIM)=' 15 ) 

CALL ASPERR 

RETURN 

END 



113 


SUBROUTINE ET PH I ( A , N A ,TT , 8 , NB , DUMMY , K DUM ) 

D IMENS I ON A ( 1 ) , B { 1 ) , DUMM Y ( 1 ) , NA ( 2 ) , NB ( 2 ) , NO ( 2 ) 
DOUBLE PRECISION A, T, TT , ANA A, TMAX, R, DUMMY 

COMMON /MAX/MAXRC 

NR=NA( 1) 

NCC=NA(2) 

NB ( 1 ) = NR 

NB<2)=NCC 
LD=NR*NCC 

IF (NR.ME.NCC.OR.NR.LT. 1 . OR . L 

NDD=2*NA( 1)*NA (1) 

I F ( KDUM .LT.NDD) GO TO 998 
NDD = NA( 1)*NA ( 1) + 1 
T =TT 

CALL NORM(A,NA,ANAA) 

T MAX= 10.01/ANAA 

K=0 

101 IF ( TMAX-T ) 103,104,104 

103 K = K + 1 

T = TT / 2**K 

IF (K-1000 ) 10 1,1 02, 102 

104 SC=T 

CALL SCALE! A,NA,A,NA, T) 

CALL UNITY(B,NB) 

11=2 

N = 35 

CALL ADD ( A,NA , B ,NB ,OUMMY ( 1),ND) 

CALL EQUATE! A, N A , DUMMY ( NOD ) ,ND) 

106 CALL MULT ( A ,NA , DUMMY (NDO) , NO , B ,NB ) 

S = 1.D0/I I 

CALL SCALE ( R »NB , DUMMY ( NDD ) ,ND, S) 

CALL ADD(DUMMY(NDD),ND, DUMMY! 1) , ND, B , MB ) 

CALL EQUATE ( B ,N B* DUMMY ( 1 ) , ND ) 

N=N-1 

IF (N) 107,107,105 

105 11=11+1 

GO TO 106 


. OR.LO.GT.MAXRC) GO TO 998 



107 IF (K) 109,108,212 

109 CALL LNCNT (1) 

WRITE (6,110) 

110 FORMAT (« ERROR IN ETPHI K IS NEGATIVE') 

112 IF ( K-l ) 213,212,212 

213 K= 1 

212 00 111 J=1 , K 
T=2*T 

CALL EQUATE IB, MB, DUMMY ( 1 ) , MO) 

CALL EQUATE (DUMMY ( 1 ) , NO, DIJMM Y ( N 00 ) , NO) 

111 CALL MULT ( 01 IMM Y ( NOD ) , NO , DUMMY ( 1 ) , NO , B , NR ) 

IXiLI 

108 CONTINUE 

S=1 .DO/SC 

CALL SCALE ( A,NA, A,NA, S) 

RETURN 

102 CALL LNCNT ( 1 ) 

WRITE (6,119) 

119 FORMAT (• ERROR IN ETPHI K=1000') 

CALL ASPERR 

RFTURN 

998 CALL LNCNT (1) 

WRITE (6,50) NA , KDUM , NDO 

50 FORMAT (' DIMENSION ERROR IN ETPHI NA=«2I6. » KDUM= ' 15 ,5X , 

1 ' KDUM ( M I M ) = ' 15) 

CALL ASPERR 

RFTURN 


END 





SUBROUTINE AUG(F,NF,G,NG,RI,NRI ,H, NH,0,NO, C,NC,Z,NZ, II ) 
DIMENSION F(1),G( 1) ,RI ( 1) ,H(1) ,Q( 1),Z (1 ) ,C( 1) 

DIMENSION NN PI (2) , NNP 21 2 ) ,NNP3( 2 ) ,NNP4(2) ,NF( 2) , MG( 2) »MR I (2 ), 

1NH(2 ),NZ(2) ,NC(2) ,NN(2 ) ,NQ(2) 

DOUBLE PRECISION F, G, RI,H,0,C,Z 

IF( (NF( 1) .NF.NFI2 ) ) .OR. (MRI (1) .NE.MRI (2 ) ) .OR. ( 

1NQ(1 ).NE.NO(2))) GO TO 995 
NNZ=2*NF( IT 

I F ( (NG( 1 ) .NE.NF ( 1 ) ) .OR. ( NG ( 2) . NE.NRI ( 1) ) ) GO TO 995 

IF ( II .EO.l) GO TO 206 

IF( (NH( 1 ).ME.NO( 1) ) .OR. ( NH( 2) .NE.NF( 2) ) ) GO TO 995 

206 TWO = 2 

N = NF(1) 

NSQ = N*M 
NZ( 1 )=MM7_ 

NZ ( 2 ) =NNZ 

NP 1 = 1 

NP2 = NP 1 + NSO 

NP3 = NP2+NSQ 

NP4 = NP3 + NSO 

CALL TRANP(G,NG,Z(NP4),NNP4) 

CALL MULT(RI,NRI,Z(NP4),NNP4,C,NC) 

CALL MULT(G,NG,C,NC ,Z (NP4) , NNP 4 ) 

IF ( II .EO. 1 ) GO TO 204 

CALL TRAMP (H.MH.Z(N°3 ) , MNP 3 ) 

CALL MULT(0,N0,H*NH,Z m , NNP 1 ) 

CALL MULTIZ (NP3 ) ,NNP3 ,7 ( 1 ) , NNP 1, Z ( NP2 ) ,NNP2 ) 

GO TO 205 

204 CALL ETHTA'T F'( 0 , NUT, Z'fNP 2 ) , 0 0") 

205 NPA I R =MOD ( N T 2 ) 

I F < NPAIR.FO.O) NPA I R = TWO 
NN ( 1 ) = N 
NN ( 2 ) = I 

GO TO ( 201 ,202 ) ,NPAIR 

201 00 300 1=1, N, 2 

NP'2 = N*( M+I-l )+l 

NTH3=TW0*N*(I-1)+N+1 







300 CALL EQUATE (Z (NIP-2 ) ,NN,7. (NTH3 ) ,NN ) 

00 302 I=2,N,2 

MP4=M»( 3»M+ I-D+l 

NTH2=TW0*N*(N+I-1)+1 

302 CALL EQUATF<Z(NP4),NN t Z(NTH2),NN) 

GOTO ( 202 , 203 ) , NP A I R 

202 DO 301 1=2 ,N,2 

M P 2 = N*(N+J-1) + 1 

NTH3 = TWD»M-!- ( I-D+N+l 

301 CALL EOIIATF ( 7. ( NP2 ) , NN ,7 ( MTH3 ) , MM) 

00 304 1=1, lv t 2 

MP4 = N* ( 3*N + I -1 ) + l 

NTH2=TW0*N*(N+I-1 ) + l 
304 CALL E0UATF(7. (MP4) , NN , 7 (NTH2),NN) 

GO TO ( 203,20 1 ) , IMPAIR 

203 no 303 1=1,0 
IJ = I+N 

no 303 J=1,N 

JJ = J + N 

L1=NMZ*( J-l)+I 

L2 = MMZ»( I vJ — 1 ) + J J 

L3=N*(J-1)+I 
7 (LI )=-F( L3) 

303 7, ( L 2 ) = F ( L 3 ) 

GO TO 1000 

995 CALL LNCMT (2) 

WRITE ( 6,50) IMF , iMG , N R I , NH , MO 

50 FORMAT (' 0 1 M ENS I ON ERROR IN A UG' , T 35 , ' M F » , 9 X , ' MG ' , 9X , ' NR I ' , 8X , 

1 * MH ' , 9 X , * NO 1 / 2 9 X , 5 ( 3X , 2 1 6 ) ) 

9 99 CALL ASPFRR 

1000 RETURN 
FMO 



117 


SUBR GUT I MB R ICAT ( PHI , NPH I , C , NO ,NC DMT , K , IMK , PT , N PT , W , K D1 IN': ) 
DIMENSION MCONT ( 3 ) , NPHI ( 2 ) ,N C ( 2 ) , WK ( 2 ) , MM ( 2 ) * MM ( 2 ) ,NPT( 2) 
DIMENSION PHI ( 1) ,C(1 ) ,K ( 1) »PT( 1) , VI 1) 

DOUBLE PRECISION PHI, C, K, PT, SUM , SUNN, A L , t 0 E T 

TWO=2 

N = MPHI (1 I/TWO 

NSO=N*N 

NSQ4=4*NS0 

NP 1 = 1 

NP2= NSO+NP1 

NP3=NSO+NP2 
NP4= NS0+MP3 

IF (kDUM.LT.NSQ4 ) GO TO 600 

IF (NPH I (2 ) .NE.NPHI ( 1) . OR.MCI 2 ) . ME . N. OR . HPT ( 1 ) .ME .M.OR.NPT (2 ) . 
1 NE.N) GO TO 600 

NPRINT=NCONT (1 ) 

NBL0CK = NC0NT(2 I/MPRINT 
M2 =NC ONT ( 3 ) 

£ REARRANGE PHI MATRIX 

NN( 1 )=N 
NN { 2 ) =1 

DO 300 1=1. N 

NTH1 =TW0#N*( I-1I+1 
NTH3= NTH1 + M 

NW1=N*( I-D + l 

NW2=NW1+N*M 

CALL EQUATE (PH I (NTH1 ) ,NN, W(NW1 ) »NM ) 

300 CALL EQUATE (PH I (NTH 3 ) ,NN,W(NW2) » NN ) 

M.M ( 1 )=TWO*N*N 
NM ( 2 ) =1 

CALL EQUATE (W ( 1 ) , NM , PH I ( 1 ) , NM ) 

DO 301 1=1, N 
NTH2=TW0*N*(N+I— 1 )+l 

NTH4= NTH2 + N 

NW3 = N* ( T W 0# N + 1 -1 ) + 1 
MW4= NW3+N*N 

r.AI I EOl IATF ( PH T ( NTH2 ) . MN « W ( NM3 ) . IMN ) 
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301 CALL EQUATE (PHI ( NTH4) ,NN,W(NW4) , NN) 

NWW=TW0*N*N+ 1 

CALL EQUATF (W (MWW ) ,NM ,PHI ( NWW ) , NM] 

C MAIN COMPUTATION 

C CALL UNITY(PT,NPT) 

DO 306 1= 1 t N 

306 K ( I) = 0. 

NT OT=0 

DO 403 1 = 1 . N R I OCK 

DO 104 J=1,NPRINT 

CALL MULT1PHI (MP3) , MPT, PT, MPT, W(l), MPT) 

CALL ADD (PHI(l), NPT , W(l), NPT, W ( 1 ) , MPT) 

CALL I N V ( W ( 1 ) , MPT, DET, W(NP2)) 

CALL MULT ( PH I ( MPA ) , MPT, PT , MPT, W(NP2), NPT) 

CALL ADD (PHI (MP2) , MPT, VMNP2), NPT, U(MP2), MPT) 

CALL MULT ( w ( NP2 ) » NPT, W(l), MPT, PT, MPT) 

SIJMNs 0. 

S U M = 0 , 

DO 202 IL=1,N 
ILP=IL+NP3 

M I L = ( I L - 1 ) * M + 1 L 

S U M = S IJ M +D A R S ( PT (MIL ) ) 

202 SUMN* SUMN+DABS(PT< NIL) -W(ILP)) 

AL = SUMM/ SUM 

DO 201 IL=1,M 
N I L = ( I L — 1 ) #M+ I L 

ILP= IL+NP3 

201 W(ILP) = PT (MIL) 

204 DO 104 M = 2 , M 

N 1 = M— 1 

DO 104 L=1,N1 
L1 = N*( L-l) + M 

I. ?=N* ( M — 1 )+). 

PT ( LI ) = ( PT( LI ) + P T ( L2 ) ) / 2 . 

PT(L2) = PT(L1) 

IF( AL-. 00001) 2Cf3 , 203 ,104 

104 CONTINUE 

MTOT=I*MPRINT 

GO TO 3 05 
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203 NT 0T=NT 0T+ J 

305 CALL MULT ( C , NC , PT , N PT , K , NK ) 

103 GO TO ( 403,400 ,401 ,402 ) , NZ 

400 CALL LMCMT (1) 

WRITE (6, 50) NT OT 
50 FORMAT ( 10 X , 1 4 , 1 4H ITERATIONS 

CALL PRNT (PT,NPT, 'P(T) • , 1) 

GO TO 403 

401 CALL LMCMT (1) 

WRITE (6, 50) NT OT 

CALL PRNT (K,NK, «K ( T) ' , 1 ) 

GO TO 403 

402 CALL LNCNT (1) 

WRITE (6, 50) NT OT 

CALL PRNT (K,NK, »K ( T) ' , 1 ) 

CALL PRNT ( PT , NPT , ' P ( T ) ' , 1 ) 

IF < AL-.00001 ) 210,210,403 

403 CONTINUE 

C REARRANGE PHI MATRIX 

210 CALL EQUATE (PHI(1),NM,W(1),NM) 

00 303 1 = 1, N 

NTH1 =TWO#N* ( I -1 ) +1 
MTH3=NTH1+N 

NW1=N»( I-l )+l 

NW2 =NW1+N*N 

CALL EQUATE (W ( NW1 ) ,NM,PH I (MTH1) , iMN ) 

303 CALL EQUATE ( W ( NW2 ) , NN , PH I ( NTH 3 ) ,N N ) 
CALL EQUATE ( PH I ( NWW ) , NM , W ( NWW ) , NM) 
DO 304 1=1, N 

NTH2=TW0*N*(N+I-1)+1 

NTH4=NTH2+M 

NW3 = N=MTW0*N+I-1 ) + l 

NW4= NW3+N*N 

CALL EQUATE ( W ( NW3 ) , NN , PH I (NTH2) , NM ) 

304 CALL EQUATE ( W ( NW4 ) , NN , PH I ( NTH4 ) ,NN ) 

RETURN 








to 

o 


600 CALL LNCNT (2) 

WRITE (6, 1600) NPH I , NIC » NP T, KDUM, NS 04 
1 600 FORMAT (» DIMENSION ERROR IN R I C AT 1 . T3 5 . ' NPH I ' . 7X , 1 NC 1 , 9X , ' N PT 1 
1 , 6 X , 'KDUM' ,3X , ' KDUM (M’l N) ' /29 X , 3 ( 3X , 2 14 ) , AX , 14, 5X , 14 ) 

CALL ASPERR 
RETURN 


END 



c 

c 


c 


c 


c 


c 

c 

c 

c 


SUBROUT INF SAM PL (PHI , NPHI , H , NH, 0 , MO , R , MR, P,NP ,K, MK , NCONT , HUM , KOIJM ) 

tttmfrtton ^"pvtit2 > , nh ( 2 rrwm ,nr ttttnp m ,nk< ts ,mcmni n irrom 


DOUBLE PRFCISION PH I ( 1 ) , H ( 1 ) , 0 ( 1 ) , R ( 1 ) y P ( I ) , K ( 1 ) , DUM ( 1 ) 

DOUBLE PRECISION SUM »S IJMN , AL 

DIMENSION OF DIJM MUST RE AT LEAST 6*N*N 

CHECK FOR CONFORMABLE MATRICES 

N=NPH I ( 1 ) 

M=NH( 1) 

NK( 1 )=N 
MK ( 2 ) =M 
NSO=N*N 
ND1 = 1 
M02 = N S0+ 1 

'KlR6=5*N50f 1 

ND3=MD2+NSO 


MS06=6»MS0 

IF (MPHI (2) .NE.N.OR.MH(2).NE.N.OR.NQ( 1 ) . NE . N .OR . MQ ( 2 ) . ME .N . OR . NR ( 1 


1) ,NE.M.OR.MR(2 ) .NE.M.OR.NP( 1 ) .NE.N.OR.NP (2 ) . NE .N. OR. KDUM .LT.NS06 ) 
2G0 TO 900 
NF I N= NCONT ( 2 ) 

NPRI NT =NC ONT ( 1 ) 


N7 = NC0NT( 3) 

START OF MAIM COMPUTATION, P { 0 ) IS INPUT DATA IN P MATRIX 

KLN=0 

J FN=0 


1 = 0 
AL = 1 . 

100 CALL MULT ( H,NH,P,NP,DUM,NZD ) 

OUM=H*P 

CALL TRAMP ( H , NH , DUM ( ND2 ) , NZD ( 3 ) ) 

0UM2= HP R I M F 

CALL MULT ( DUM , NZD , DUM ( NO 2 ) ,N ZD( 3 ) , DUM ( ND3 ) , NZ D ( 5 ) ) 
DUM3=H*P«HPRIMF 

CALL ADD ( DUM ( ND3 ) , N Z 0( 5 ) , R , NR , DUM, NZ 0 ) 

DUM=H*P*HPR IME+R 

CALL PSEUDO ( DUM, NZD , DUM ( ND 2 ) , N Z D ( 3 ) , DUM ( ND3 ) , KDUM-3*NS0 ) 

DUM2= INVERSE 

CALL TRAMP ( H , NH , DUM , N Z D) 

DIJM=HPRIMF 

CALL MULT ( DUM , N ZD , DUM ( ND 2 ) , N Z.D( 3 ) , DUM ( NQ3 ) , NZ D ( 5 ) ) 



£ CALL MULT ( P , NP , DUM ( ND 3 ) , NZ D ( 5 ) , DUM , NZD ) 

^ C niJM = A 

110 CALL MULT ( PH 1 1 MPH 1 1 PI IM , NZ Dt K> MK ) 

CALL MULT ( niJM , NZO ♦ H , NH ,DUM ( ND2 ) , NZD ( 3 ) ) 

C DUM2=A*H 

CALL MULT ( HUM ( ND 2 ) , N ? D( 3 ) , P »N p , HUM , NZ 0 ) 

C DUM=A*H*P 

CALL SURT (P,MP,DUM,N7.D t OUM,NZD> 

_C nUM= P-A*H*P 

CALL TRAMP ( PHI ,MPHI , DUM ( ND 2 ) ,NZD ( 3 ) ) 

CALL MULT ( DUM , NZD , DUM ( ND2 ) , NZ D ( 3 ) , DUM( MD3 ) , NZD ( 5 ) ) 
CALL MULT ( PH I ,NPHI , DUM ( ND3 ) » NZD ( 5 ) , DUM , MZ D) 

T UUM= T R1 (F-A*H*P )MHlPKlMb 

CALL ADD (niJM,NZD»0,M0 t P,NP) 

00 120 M=2,M 
Ml = M-l 

00 120 L = 1 1 N 1 

L1 = N*( L-D + M 

L 2 =N* ( M- 1 ) +L 

P ( LI ) = ( P (LI )+P (L2 ) ) /2. 

120 P (L2) = P (LI) 

130 IF (I .EO.O) GO TO 150 
SUM=0. 

Sl)MN=0 . 

DO 140 I L = 1 » N 
IJ= ( I L — 1 ) *N+ I L 

SUM=SUM + OARS(P( I J) ) 

NDL =MD6 + I L 

140 SUMN= S U M M + 0 A R S ( P ( I J ) —DIM ( NDL ) ) 

AL = SUM N/ SUM 

150 DO 151 IL=1,M 

I J= ( I L — 1 ) *N+ I L 

NDL= ND6+ I L 

151 DUM (NDL) = P ( I J ) 

I L ST= I 

1 = 1 + 1 

IF (AL.LE.. 00001) GO TO 300 
IF ( I .GE.NFIN) GO TO 310 
C INTFRMEDIATF print 



IF (I .LT.MPRIMT ) GO TO 100 
MP RI MT = N PR I MT + NC ON T ( 1 ) 

- 15 ? . fiQ . TO ( 17 0.156.1 55.1551.Ni7 

155 CALL LNCNT ( ? ) 

WRITE (A, 1152) ILST 

115? FORMAT ( * 0 S T F P NUNBFR= 1 « I 4. 1 IN S A M P l_ 1 ) 

CALL PR NT (K,NK,4HK ( T) ,1) 

GOTO (170,170,170, 156), MZ 

156 CALL L NCNT ( 2 ) 

WRITE (6,1152) I 

160 CALL PR NT ( P , NP , 4HP ( T) , 1 ) 

170 IF (JFN.EQ.O) GO TO 100 

RETURN 

300 CALL LNCNT ( 2 ) 

WRITE (6,1300) 

1300 F OR M AT ( ' OP NO L ONGER CHANGING, LX I I FROM SA^PL ' ) 

310 JFN=1 

GO TO 152 

" OITCrTTAlL ' LMLN'I (2) 

WRITE (6,1000) MPHI ,MH,NO,NR,NP ,K0UM,MS06 
1000 FORMAT ( ' 0 T M E N S I OM ERROR IN SAMPL',T35 , 3 X , 4HN P H I , RX , 2HMH , 

19X ,2HNQ » 9X , 2HNR , 9X » 2HNP , 5X » 4HK0UM, 3X » 9H KOUM ( M IN ) /2 9X , 5( 3X , 2 14 ) , 
23X, 14, 5X, 14) 

CALL ASPERR 

RETURN 

END 


N> 

U> 



SUB ROUT INE TRNS I ( F , MF , G , NG , J , NJ , R , NR t K , NK , H , NH , X , NX , T , DUMMY , KOUM ) 


(O 

4s- 


c 


c 


DOUBLE PRECISION F ,G , J, K ,H , X ,T , Y , R , DUMMY 

DIMENSI ON F( 1) ,NF( 1 ) ,G ( I) ,NG( 1) , J( 1 ) ,NJ ( 1 ) ,K( 1 ) ,NK ( 1 ) ,H( 1 ) ,NH( 1 ) 

1X( 1) ,NX( 1 ) ,T( 1 ),R ( 1) ,NR( 1 ) t DUMMY ( 1 ) 

DIMENSI ON NNF ( 2 ) , NNG(2 ) , NU ( 2 ) , NN X ( 2 ) , NY ( 2 ) 

DATA STAR/'*'/ 

I F ( NF ( 2 ) • NE • NG ( 1) .OR. NJ ( 2 ) . NE . NR ( 1 ) .OR. NK ( 2 ) .NE.NX ( 1 ) .OR. 

IITJIl J.Nt.NK "("II '.UR. NRTTT.Nb. NX ( 2T — .'UK'. NH( 2 ) .NE .NX ( I ) .UIH 

2NF(2) .NE.NX(l)) GO TO 999 
MAX = NF ( 1 ) * ( NF ( 2 ) +NG(2)+ 1) +NH(1)+NK(1) 

IF ( KDUM .LT. MAX) GO TO 910 
II = 1 

NSO =N F ( I 1 ) *NF (II) 

NX4 = NSO *4 

I F ( KDUM .LT. NX4 ) GO TO 900 

TRNS I PROGRAM 
N2 = II + NSO 

N3 = N2 + NSO 

N4 = N3 + NSO 

L 3 = N2 +NF( I1)*NG(2) 

L4 = L 3 + N J ( I 1 )*NR( 2 ) 

L5 = L4 + NH ( 1 ) 

L6 = L5 + NJ( I 1 )*NR( 2 ) 

T 1 =T( 1 ) r 

N = ( T ( 2 ) + • 5*T 1 ) / T 1 

NXR =NX (ID 
LAST = L 6 -II 
TT = T ( 4 ) 

CALL PRNT ( F » NF » ' F ',1) 

CALL EAT(F t NF,Tl,DUMMY ( N3 ) , NNF, DUMMY ( N4 ) , NNF, DUMMY (II ),KDIJM) 

100 FORMAT ( 1 HO 1PSE16.7) 

CALL PRNT ( DUMMY ( N3 ) * NF, 'EAT ', 1) 

CALL EQUATE ( DUMMY ( N3 ) , NF , DUMMY ( II), NNF) 

CALL P RN T ( DUMM Y ( N4 ) , NF, "'INT ', 1) 

CALL MULT (DUMMY ( M4) , NF , G , NG , DUMMY ( N2 ), NNG ) 

CALL MULT ( J,NJ,R,NR,DUMMY( L3) ,NU) 

CALL LNCNT (100) 

CALL LNCNT (3) 
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WK1 , b|6t bU , 

50 FORMAT ( 1H0 'TRANSIENT RESPONSE, * INDICATES CONTROL CHANGES') 

WR I TE ( 6 , 51) NXR , NH ( I 1 ) , NK ( 1 1 ) 

51 FORMAT! 1H0 AX, 'TIME FIRST', 13,' ELEMENTS CONTAIN X, NEXT', 13,' 
1 ELEMENTS CONTAIN Y = HX, LAST', 13,' ELEMENTS CONTAIN U =JR -KX' 
2) 

201 N0=1 

IU = 11 

CALL MULT ( K, NK, X,NX, DUMMY (L5), NU) 

CALL SUBT (DUMMY! L3) , NU, DUMMY ( L b ) , NU, DtJMM Y ( L 5 ) , NU ) 

203 CALL MULT ( H,NH, X ,NX , DUMMY ( L A ),NY) 

IF (IU .NE. ID GO TO 205 

WR I T E ( 6 , 101) STAR,TT, (X! IP), I P = 1 , NXR ) , ( DUMM Y ( I P ) , I P =L A , L AS T ) 

101 FORMAT ( 1 HO A1 ,F8 .2 , 1P7D1 5 .7/ ( 1 OX , 1P7D15. 7 ) ) 

GO TO 206 

205 WR ITE ( 6, 102) TT, (X ( IP) , IP=1,NXR) , (DUMMY! I P ) , I P = L A , L AST ) 

102 FORMAT ( 1H0 IX , F 8 . 2 , 1P7 D15 . 7/ ( 10X , 1 P 7D1 5 . 7 ) ) 

206 IU = II + IU 

IF (TT .GE. DABS( T ( 3 ) ) ) RETURN 

TT * TT + T1 

202 CALL MULT ( DUMMY (II), NNF, X, NX, DUMMY(L6 ),NNX) 

CALL MULT (DUMMY ( N2 ) ,NG» DUMMY ( L 5 ) , NU, X, NNX ) 

CALL ADD(X»NX, DUMMY ( L6 ),NNX, X,NNX) 

I F ( NO ,GE . N) GO TO 201 

NO = NO + 1 
GO TO 203 


C DIMENSION ERROR DIAGNOSTIC 
999 WR I TE ( 6 , 990) 

990 FORMAT ( 1H0 'DIMENSION ERROR IN TRNSI ' /25X, ' COL SIZE OF 1ST MATRIX 
1 ROW SIZE OF 2ND MATRIX') 

WR ITE ( 6 » 99 1 ) NF ( 2 ) , NG(1) 

WR I TE ( 6 ,992 ) NJ(2),NR( 1) 

WR ITE ( 6 , 993 ) NK(2),NX(1) 

WRITE(6,995) NH(2),NX( 1) 

WR IT E ( 6 ♦ 996 ) NF(2), NX(1) 

WR I TE ( 6 »99A ) N J ( 1 ) , NK ( 1 ) , NR ( 2 ) , NX ( 2 ) 



126 


991 F ORMAT ( 1HO 'INTEGRAL (EXP(F.T)) G ' 113, 20X, 18) 

992 FORMAT ( 1 HO 'J R' 17X , 1 15 , 20X , I 8 ) 

993 FORMAT! 1HO »K X* 17X , I 15, 20X ,18) 

994 FORMAT ( 1 HO 'JR -KX ROW SIZE OF J IS'I5,3X,'OF K IS ' , I 5 ,3X , 'COL 

1 SIZE OF R IS' 15, 3X, 'OF X IS' 15) 

995 FORMAT ( 1 HO 'H X' 17X , 1 15 , 20X ,18) 

996 FORMAT ( 1HO 'EXP(F.T) X* 10X , I 15, 20X , 1 8 ) 

GO TO 1000 

900 WR I fE ( 6 , 52) nx4,k0um 

52 FORMAT { 1 HO 'DUMMY MUST BE DIMENSIONED AT LEAST' 14, • X 1' 'BUT IS 
1 DIMENSIONED ONLY' 14,' X 1') 

GO TO 1000 

910 WR I TE ( 6 , 52) MAX , KDUM 
1000 CALL ASPERR 
RETURN 

END 
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SUBROUTINE PSEUDOt A, NA, B, MB , HUM , KDIJH ) 

DIMENSION A ( 1 ) , B( 1 ) ,DUM (1) ,NA ( 2) ,NR ( 2 ) , I p ( 4 ) , DE P ( 3 ) 

DOUBLE PRECISION A , R , DIJM , DE P 

NNW = 3*NA( 1 ) *NA (2) 

DO 10 1=1,2 
DEP ( I ) = 0 . 0 

10 IP(I)=0 

20 I P ( 3 ) =NA ( 1 ) 

IP( 4 ) =N A ( 2 ) 

NB( 1 )=NA(2) 

NB( 2 ) =NA ( 1 ) 

IF ( KDUM.LT.NNW) GO TO 999 

NEE = NA( 1 )»NA(2) + 1 

ND=2*NEE - 1 

C IF A IS (1,1) MATRIX ROUTINE INVERTS All) AMO PUTS IT IN R(l) 

LE (NAtl) .FQ. 1 .AND . N A.L2-) .FQ.-1) GO T O B OO 

CALL PSEtMA , R , DUM , DUM ( NEE ) , DEP , I P, DUM ( NO ) ) 

GO TO 1000 

600 B( 1 ) = 1 . /A ( 1 ) 

GO TO 1000 

999 WRITE(6,500) KDUM,NNW 

500 FORMAT (» DIMENSION ERROR IN PSEUDO K OUM= » I 5 , 3X , 1 K DUM ( M IN ) = 1 I 5 ) 

1000 RETURN 
END 



C SUBROUTINE DFCGEM DECOMPOSE A Rb A L G F N E R A L MATRIX 

SIJBR OUT IMF O EGG EM ( R , NR ,0,M0,,H, NH , HUM , KDIIM ) 

D I MENS I ON NR (2 ) ,NH(2 ) ,NG( 2 ) ,NO( 2 ) , NA ( 2 ) ,KP ( 4) 

COMMON /MA X/ Pi AX PC 

DOURLE PRECISION R ( 1 ) , G ( 1 ) , H ( 1 ) , RUM ( 1 ) , OCM ( 3 ) 

NJ=MR( I ) *NR ( 2 ) + 1 

CALL TRAMP ( R , NR , DIJM ( M J ) , ND ) 

I F ( MR ( 1 ) .GT.NR ( 2) )GO TO ID 

CALL MULT ( R » NR » DUM ( N J ) , NO »f)UM ,M A ) 

ICS=1 

DO TO 30 

10 CALL MULT ( DIJM ( M J ) , MD , R , MR r DUM, NA ) 

ICS = 2 

GO TO 30 

C ENTRY DECSYM DECOMPOSE A REAL SYMETRIC MATRIX 

FMTRY DFCSYM ( R , NR , G . MG , H , MH , OHM , K DUM ) 

ICS=0 

IF ( MR ( 1 ) .ME. NR (2 ) ) GO TO 601 

CALL EQUATE ( R , MR , DUM , M A ) 

C DUMMY STORAGE MUST RE AT LEAST 7*NA(1)**2 

30 M = N A ( 1) 

MM 7= 7 =:-M»M 

IF ' (KDUM.LT.MM7 ) GO TO 601 
KP( 1 ) =0 

KP ( 2 ) =0 

K P ( 3 ) = M 
KP ( 4) =M 

P CM ( 1 )=0. 

DCM (2 ) =0 . 

MS=M*M 

N2 = MS+1 

N3=N2+MS 

M4=M3+MS 

N5=N4 + MS 

N 6 = N 5 + M S 
M7=N6+MS 

CALL DEC DM (DUM.DUM(N7 ) , D IM ( N4 ) < DUM ( N5 ) . niJM ( N6) t OCMt KP« DUM( M2 ) ) 
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CALL TRAMP ( DUM ( N3 ) , NA , DUM ( N6 ) « ND ) 

IF IICS.E0.2) GO TO 100 

CAM MJLLL T (DllMIN?) .NA.DUM(N6).N0.H.NH) 

NH(2 )=KP( ?.) 

IF ( ICS.FO. 1) GO TO 60 

CALL TRAMP (EUNH,G,NG) 

GO TO 150 

60 CALL MULT ( DUM ( M5 ) , N A ,H , NH ,6 ,NG ) 

CALL TRANP ( G ,NG , DUM ( N6 ) , NO ) 

CALL MULT ( DUM ( N6 ) , ND , R , NR , G , MG ) 

GO TO 150 

100 CALL MULT ( DUM ( N2 ) » M A , DU M( N6 ) , MU , H , MH ) 

NG( 2 ) =KP ( 2 ) 

CALL TRANP (H T NH,G,MG) 

CALL MULT (G>NG,DUM ( N5 ) »NA t H,NH ) 

CALL TRANP (H,NH ,DUM { N6 ) » NO ) 

CALL MULT ( R ,NR , DUM ( N6 ) , NO, H , NH ) 

150 DUM (N6)=KP( 2 ) 

RETURN 

601 CALL LNCNT (1) 

WRITE ( 6_f 1601 ) NR»KDUM,MM7 

1601 FORMAT ( • DIMENSION ERROR IN DECSYM ( DECGFN ) NR='2I6,3X, 

1 'KDUM='I4,3X, 'KDUMIM IN)=' 14) 

CALL ASPERR 

RETURN 

END 



w 

o 


SUBROUTINE RE AD 1 ( A, NA, NZ ,NAM) 

COMMON /MAX/MAXRC 
DIMENSION A ( 1 ) ,N A ( 2 ) , NZ ( 2 ) 

DOUBLE PRECISION A 

IF (NZ(l).FO.O) GO TO 410 
NR = NZ( 1) 

NC=NZ ( 2 ) 

NLST=NR*NC 

IFIMLST .GT. MAXRC .OR. NLST .LT. 1. OR.NR.LT. 1) GO TO 16 

DO 400 I s 1, NR 

400 READ (5,101) ( A ( J), J = I, NLST, NR) 

NA ( 1 ) =MR 

NA( 2 ) =NC 

410 CALL PRMT ( A , NA , NAM , 1 ) 

101 FORMAT (8F10.2) 

RETURN 

16 CALL LNCNT(l) 

WRITE (6,916) NAM,NR,NC 

916 FORMAT (' FRROR IN REA01 MATRIX 'A4,' HAS NA=',2I6) 

CALL ASPERR 

RETURN 

END 



SUBROUTINE ASPERR 

DATA I /10/ 

CALL TRACE 

C ERRTRA IS THE 360/67 TRACE ROUTINE TRACE IS FOR TSS 

C CALL ERRTRA 

C THIS IS AN INSTALLATION DEPENDENT SUBROUTINE 

C SUBROUTINE ERRTRA IS A SUBROUTINE SUPPLIED BY THE AMES OPERATING 

C SYSTEM TO PROVIDE AN ERROR WALKBACK 

C THE STATEMENT CALL ERRTRA SHOOLD BE EITHER 

C 1) CHANGED TO MATCH THE USERS OPERATING SYSTEM, 

C OR 2) DELETED ALTOGETHER. 

1 = 1-1 

IF (I.GT.O) RETURN 

1 = 10 

WRITE (6,100) 

100 FORMAT C TOO MANY ERRORS. EXIT CALLED 1 ) 

CALL EXIT 

RETURN 

END 
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BLOCK ..QAT A 

COMM ON /FORM /NE PR, FMT1 ( 6) , FMt2 (6) 

COMMON/ LINE S/NLP, LIN, TITLE ( 23 ) 

COMMON /MAX/MAXRC 

DATA MAXRC/6400/ 

C- MOTE NLP MO. LINES/PAGE VARIES WITH THE INSTALLATION. 
DATA LIN, ML P/1, 45/ 

DTTA NEPR ,FMT1 / 7, M 1P7D16.7) ' / 

DATA FMT2/' (3X,1P7D16.7) '/ 

DATA TITLE /19»' »,*VASP PROGRAM '/ 


END 
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C- SUBR TO COMPUTE PSEUDO-INVERSE OF GENERAL MATRIX, RETURN FINAL PIVOT 
C... NOTE IMPLIT STATEMENTS MUST BE -FIRST- CAM RE REPLACED BY TYPE 

IMPLICIT REAL#8 (D), INTEGERS (0) 

COMMON /MAX/MAXRC 

C DOUBLE PRECISION IS THE ONLY THING ESSENTIAL. 

INT EGER#2 M 

DOUBLE PRECISION A,B,C,EE, D 

DIMENSION A(400),B(400),C( 400) , EE ( 400 ) , 0(2000), 

1 KRV(4), 

2 DE P( 3) , DPR ( 2 ) , IP(4), JP(5) 

DATA ICC, DFZO /Z40000000, O.DO / 

EQUIVALENCE (DD I ,FDI , IDD) , ( DMX , FMX ) 

EQUIVALENCES I ,DSUM) , ( DFZ 0, FZRO , I Z , OZ ), ( 0LL,0R1 ) , ( KRV ( 1 ) , KRC ) , 

1 (KRV(2)tKRC2),(KRVO),KRC3),(KRV(4),KRC4) 

OPS = 1 
GO TO 1000 

FNTRY PSEU PtA.B.C.EE.DEP. IP, D) 

QPS = IZ 


I P( 4 ) = IPO) 

1000 CONTINUE 

DPI = DEP(l) 

EF2 = SNGL ( DEP (2 ) ) 

C- SET DEFAULT VALUES OF TOLERANCES 

IF(DEP(1) .EQ. DFZO) DPI = 2.D-6 

IFIEF2 .EO. FZRO) EF2 = 1.0 
NCA = I P ( 4 ) 

C NUMBER OF ROWS OF ORIGNAL INPUT MATRIX 
QR = IPO) 

C- SET SW FOR =0, DO ALL STEPS, N0T=0, THEN WANT RANK ONLY. 

ONT = QR*NCA 

C- TEST DIMENSIONS INPUT FOR REASONABLENESS. 

I F ( QNT .LT. 2 .OR. QNT .GT. M AXRC .OR .OR . LT . 1 ) GO TO 691 
C- IF DIMENSIONS ABSURD, PSEU ERR EXIT 1. 

QDCM = I P ( 1 ) 

QITR = QDCM 
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IF ( QDCM .LT. QZ) QITR = QDCM +1 
NR = QR 
OIT = IZ 



Z c- TEST TO SEE IF SYMMETRI ZATI ON IS NEEDED. 

* I F ( OPS ) 16, 150, 16 

C- TEST TO FIND SMALLER DIMENSION OF MATRIX. 

16 I F ( OR - MCA) 18,18,19 

19 MR = NC A 

OX = 1 

0R2 = OR 
OLL = OR 

OTP = IZ 

GO TO 170 
18 CONTINUE 

OX = NR 
0R2 = 1 
OLL = NCA 
OTP = 1 
170 CONTINUE 

C- SET ROW-COLUMN LIMIT TO APPROPRIATE CASE, EITHER ROW OR COLM DIMENS. 

DO 181 I =1. NR , 

DO 181 K =1, NR 
DSUM = DF7.0 

DO 183 J = 1. OLL 

ON = OX*J - OR 
L = ON + 0 R 2 * I 

M = ON + 0R2*K 

183 DSUM = DSUM + A(L)*A(M) 

C- SUM OF A( I , LL ) * A ( K, LL ) , , LL RUNS 1 TO ROW OR TO COLM LIMIT 

C . R = A»A TRANS HAS COLM LIMIT, B = A TRANS * A HAS ROW LIMIT. 

L = ( K-l) *NR + I 
181 R ( L) = DSUM 

GO TO 188 

C- HERE MOVE A TO R. A IS ALREADY POSITIVE DEFINITE. 

150 DO 151 L = 1, ONT 

151 B ( L) = A ( L ) ■ 

C-. FORCE SYMMETR I ZA T I ON OF R, TO COMPENSATE FOR ROUND-OFF, MULT IPL IC . 
18.8 DO 189 I =1, NR 

DO 189 K =1, NR 

~t BTTTkT = B l K , I ) = 1/2 ( B ( I , K ) + B ( K, I > ) 

L = I + (K-1)*NR 
M = K + ( 1-1 )*NR 
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DSUM = ( B ( L ) + B(M) ) * 0.5D0 
B ( L) = DSUM 

_L£2 RIM) = nstiM 

C HERE SET UP CALL -INITIAL- OF ANDRA. ONLY COMES HERE ONCE PER MATRIX. 
ONT = NR*MR 

KRC = ONT 

KRC2 = ONT + KRC 

KRC3 = ONT + KRC2 

KRC4 = ONT + KRC 3 

C-* OMIT SAVING OF B, IF RANK ONLY AND NO ITERATION 
IF ( IP( 2) ,NE. IZ .AND. OITR .EP. IZ) GO TO 200 
DO 1891 I =1, ONT 
1891 D C I ) = B(I) 

200 CONTINUE 

C- SEARCH DIAGONAL OF INPUT FOR LARGEST ELEMENT. USE TO DEFINE EL. PT . 

0R1 = NR + 1 
L = 1 

DMX = DFZO 
M = IZ 

DO 23 I = 1, NR 

DPI = DABS ( B ( L ) ) 

IFIFMX .GE. FDI ) GO TO 23 
M = L 

FMX = FDI 

23 L = 0R1 + L 

I F( M .EQ. IZ) GO TO 692 

C- SET T OL FRANC F FOR ANDRA LIMIT OF SIZE OF DIAGONAL. 

C TOLERANCE OF ZERO IN ANDRA CALL. 

DPR (1) = DABS ( DPI * B ( M ) ) 

C ~ ASK FOR ALL ROWSt DONE IN 1 CALL 

JP ( 1 ) = IZ 

C- JP2 FIRST TIME INITIALIZATION FOR ANDRA 

JP ( 2 ) = IZ 

I F ( 0 1 T .NE. OZ) GO TO 561 
JP{4) = NR 

M = IZ 

SOLD = - EF2 

C — HAVE FINISHED PRELIM. PART 

C INITIALIZATION FOR ANDRA ( DI AGONA LI Z AT ION ) NOW COMPLETED. 
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CALL ANDRA -TO DIAGONALIZE SYMMETRIC MATRIX. 

C CALL ANDRA REDUCES ROWS BY MODIFIED GAUSS METHOD, USING SORT(PIVOT). 

30 CONTINUE 

TFTTTT IK .E ' CJ. "OZ) GU I U 3 7 

C- SAVE OLD VALUES IN CASE PIVOT IS REJECTED, UNDER ITERATION OPT. 

DO 31 L =1, ONT 

J = KRC + L 
K = KRC2 + L 

D< J> = B (L ) 

31 D ( K ) = C(L) 

32 CALL ANDRA (B,C, DPR, JP) 

JP ( 1 ) = QITR 

IR = JP ( 3 ) 

C- CHECK COMPLETION- IS MATRIX ALL D ONE IS MATRIX INVERTIBLE.. 

I FI QITR .EQ. IZ .OR. IR .EQ. NR .AND. Q IT ,E0. IZ) GO TO 700 

CHECK IF ITERATING WITH RHO TEST OR NOT 
C* OUT IF NO ITERATION OR NO NEW PIVUT FOUND 

C- OMIT ITERATION CALCS. IF NO NEW PIVOT. DECREASE TOLERANCE 

I F ( J P ( 5 ) .FO. M) GO T 0 220 

C COMPUTE RHO FOR ESTIMATING THRESHHOLD TO STOP SS IS RHO 

SS = (BDNRM(NR,C,EE,Q,KRV) + BD NRM ( -NR , C , EE , D , KR V ) ) *EF2 ** IR 

C WHY ONLY SNGLE PREC./THIS IS ONLY A ROUGH TEST TO STOP ITERATION. 

C THAT-S WHY. SIMILARLY, OTHER USES OF SINGLE PREC. 

IFISOLD .LT. SS .AND. SOLD .GT. FZRO) GO TO 650 

C- IF SUBSTANTIAL IMPROVEMENT TRY AGAIN, 

C, OTHERWISE QUIT, RETURN THE A PSEUDO INVERSE, EVEN IF OFF. 

220 CONTINUF 

OIT = on +1 
SOLD = SS 

C/ SAVE PREVIOUS ROW IN WHICH A PIVOT WAS FOUND 

M = J P ( 5 ) 

I F ( OIT .EQ. NR) GO TO 700 

C- PUT IN SMALLER TOLERANCE IN CASE DIAGONAL TOO SMALL OTHERWISE. 

DPR (1 ) = DPR (2 ) * 2.0-5 
C- TRY TO REDUCF 1 MORE ROW. 

I F ( I R ~ NR) 30 , 700, 606 

650 CONTINUF 

C* RESTORE B AMD C TO THEIR PREVIOUS VALUES. THE LAST PIVOT HAS BEEN 
CREJECTED (BACK-TRACK), WHILE ITERATING. 
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JP( 3 ) = JP(3) -1 

DO 653 I =1, ONT 

J = KRC + I 

K = KRC2 + I 
B ( I ) = D(J) 

653 C ( I) = D(K) 

700 CONTINUE 
IR = JP ( 3 ) 

M = IZ 

C- HERE WISH TO REPLACE MARKERS IN DIAGONAL WITH LEGITIMATE l.DO 
L = 1 

DO 7 04 I = 1, NR 

DDI = B(L) 

I F ( IDD) 70 1, 702, 701 

701 IF ( IDD .ME. ICC) GO TO 7101 

B ( L ) = l.DO 

GO TO 704 

C AT 7101 FORCE SMALL TRASH TO ZERO. 

7101 B ( L ) = DFZO 

702 M = 1 

704 L = 0R1 + L 

C - IF ALREADY TRIED ANOTHER REDUCTION, TO GET MATRIX IN -UPPER- UIAG. 
COR OMIT PART OF CALCULATIONS IF ONLY RANK IS DESIRED. 

IF( I P ( 2 ) .NE. IZ) GO TO 877 
M OUCH ^URRrESSES La St PHASE IF DECOM WAS CALLtR.. 

IF(M .LT. 1 .OR. ODCM .LT. OZ) GO TO 80 

C BELOW HAVE SING. MATRIX THAT NEEDS FURTHER WORK. 

C- HAVE MATRIX DIAGONALIZED WITH IS, OS INTERSPERSED (A IS SINGULAR) 

C- RE-00 TO GET PS-INV THAT MOVES ALL IS OF DIAGONAL TO UPPER LEFT DIAG. 
C .TO COMPUTE U MATRX AS IN ASP, FOR TRANSFORMING ORIG B IN SINGULAR CASE 

rirr 

DO 527 I =1,NR 

DO 525 J =1, NR 

K = (J-1)*NR + I 

IF ( B( 4 ) ) 521,522,521 

522 C ( K) = ~ C ( K ) 

C (L ) = DFZO 

GO TO 525 

521 C ( K ) = DFZO 



CIL) - l.DO 
525 CONTINUE 

52 I L = OR1 + L 

SAVE RANK SO FAR, SHOULD BE SAM E SIZE AFTER RE- INVERSION 
0R2 = IR 

DO 54 I = 1, NR 

DO 54 K =1, NR 
DSUM = DFZO 

ON = <K-1)»NR 

DO 53 J =1, NR 
M = ( 1-1 >*NR + J 

L a ON + J 

53 DSUM * C(M)*D(L) + DSUM 
L * ON + I 

EE( L ) = DSUM 

54 CONTINUE 

DO 56 I =1, NR 

DO 56 K =1, NR 

DSUM = DFZO 
ON = { K-1)*NR 

DP 55 J =1 . NR 

L = (J-1)*NR + I 
M = ON + J 

55 DSUM * EE ( L ) ^C ( M ) + DSUM 
L a ON + I 

B (L ) = DSUM 

56 CONTINUE 

C, SET UP FOR SECONDARY ANDRA CALL NO ITERATION J-P4 = MR 
0 I T =1 

C GO FIND LARGEST DIAG. ELEMENT AGAIN 

GO TO 200 
561 JP( 3 ) = IZ 

CALL ANDRA ( BtEE ,PPR, JP ) 

IR = J P ( 3 ) 

C- TEST FOR A CHANGE IN RANK ... ERROR 

I F ( QR2 ~ IR) 693, 56ft, 694 

568 CALL T TRM ( NR , EE » D) 

C- TRANSFORM C SHARP IN D.. BS = ( ( U ) * D *(U TRP ) ) 

DO 58 I =1, NR 


DO 58 K =1, NR 
DSUM = DFZO 

ON = (K-1)*NR 

DO 57 J =1, NR 
M = (J-1)*MR + I 

L = ON + J 

57 DSUM = C(M)*D(L> + DSUM 
L = ON + I 

B ( L ) = DSUM 

58 CONTINUE 

DO 60 I = 1 » MR 

DO 60 K = 1. MR 

DSUM = DFZO 
DO 59 J =1, NR 

ON = (J-1)*NR 

M = ON + K 
L = ON + I 

59 DSUM = B ( L ) »C ( M ) + DSUM 

L = (K -1)*NR + I 

E E ( L ) = DSUM 

60 CONTINUE 

C- MOW RE-ENTER MAIM SEQUENCE WITH PS-INV. IN FE. 

GO TO 808 

C GO FIX UP B PS U EDO- INVERSE. PRESUMABLY HAVE n I AGONAL I ZED 

C HAVE DIAGONALIZED WITH ALL IS IN UPPER LEFT 

C- HERE WE HAVE FINISHED DIAGONALIZ. WANT TO GET PSUEDO INV. IN B. 

870 I F ( QDCM .LT. OZ) GO TO 877 

C NEED TO SAVE DIAGONALIZED B FOR USE BY DECOM CALL (ODCM MEG. FLAG) 

DO 871 I = 1 , ONT 

C- A WAS SYMMETRIC. JUST MOVE EE TO B RETURN FROM PSEUP ENTRY 

871 B ( I ) = EE ( I ) 

GO TO 877 

80 CONTINUE 

C MOW FORM (T TRP ) * T = APPROX B SHRP PSUEDINV IN MAT RX EE 
CALL T TR M(NR,C,EE) 

808 I F ( QPS .EO. OZ) GO TO B70 

IF(OTP) 819, 819,818 

h- C HERE B = (A TRANS)*E = (A TRP)*(A*A TRPJ-SHRP NRA .LE. NC A 

818 DO 8181 I = 1 , NCA 
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DO 8181 J = 1,NR 
DSUM = DFZO 

ON = (J - 1 ) *NR 

DO 8182 K = 1,NR 
L = ( 1 -1 )*OR + K 

H = ON + K 

8182 DSUM = DSUM + A(L)*EE(M) 

L = (J-1)*NCA + I 

8181 B ( L) = DSUM 

GO TO 877 

819 DO 8191 I = 1 »NR 

DO 8191 J =1. OR 

DSUM = DFZO 

DO 8192 K = 1,NR j 

L = (K -1)*QR + J 

M = (K -1)*NCA + I 
8192 DSUM = DSUM + A(L)*EE(M) 

C- NOTE NCA IS USED* BECAUSE A-SHARP IS TRANSPOSED IN DIMENSIONS 

L = (J -1)*NCA + I 

C HERE B = EE (A TRANS) = (A TRP*A)-SHRP * (A TRAMS) NR A .GT. MCA 

8191 B ( L) = DSUM 

C- HERE GET READY TO RETURN 
877 CONTINUE 

C- MOVE RANK TO RETURN PARAMETER 

I P( 2 ) = IR 
DEP (3) = DPR ( 2 ) 

C. ABOVF RETURN FINAL PIVOT FROM ANDRA ALG. PI AGONAL I ZAT ION 

RETURN 

691 CALL LNCNT(l) 

WRITE (6,1691) OR, NCA 

1691 FORMAT <• DIMENSION ERROR IN PSEU NA='2I6) 

GO TO 1700 

692 CALL LNCNT(l) 

WRITE ( 6,1692) 

1692 FORMAT (' ERROR IN PSEU - DIAGONAL ELEMENTS OF MATRIX=0') 

GO TO 1700 

693 CALL LNCNT(l) 

WRITE (6,1693) 

1693 FORMAT <« ERROR IN PSEU RANK HAS DECREASED COMPUTATION ENDED*) 


iii ii iii ii mu in i 



— O 
— >fr cc 
hct q: r 

■*— vC LLl Qd 

O j— • — 1 a: ► 

02 - - uj oo; 

o o ^ cl sot 

*-+ Z. *s> in- 

I h- <1 

O LLI C O 


a: 

— o 

sO cc. 

O QC 
vO LU 

i— I QC 

► — a: 

UJ 

— * CL 


JM J _J hhOC -J J— Q 
OcqcOccDc qc O < UJ 2 

C u2 U. U OO 2 1L u C^IU 


141 



S FUNCTION BDNRM(NR»CT«EE«D,KRV) 

I NTEGER*2 OF 

DOUBLE PRECISION CT,EE,0, AN,BR, DFZO, DSUM 

DIMENSION CT(40Q),EE(400), 0(2000), NV(2), KRV(4) 

C- D HOLDS 5 MATRICES. THE FIRST AMO THE LAST 2 ARE USED HERE 
DIMENSION PPP ( 2 ) 

EQUIVALENCE ( AN , FN ) , (BR,FR) 

DATA DFZO /O.DO/ 

Ct EQUIVALENCES BELOW JUST TO SAVE STORAGE 

EQUIVALENCE (DFZO, I Z ), (AN,DSUM) , ( BR,PPP ( 1 ) ,1 ) , ( PPP ( 2 ) ,K) , 

1 (NV(1) ,L), (NV(2) ,M), ( IR,NL) 

C TEST,, IF NR NEG . , THEN TRANSPOSE ROLES OF U AND (CT TRAMS *CT) 

OF = NR 
Kbi = KRV ( 3 ) 

KD4 = KRV ( 4 ) 

IF(NR) 10.10. 20 

ENTRY TTRM ( NR »CT, EE ) 

C TO DO T TR 4 T ONLY ENTRY TTRM 

QF = I Z 

GO TO 20 

10 NR = -NR 

20 IR = NR 

DO 30 I = 1, IR 
LL = ( 1-1 ) * I R 

PQ...1Q. . K = . 1 1 J R 

DSUM = DFZO 
KK = (K -1 ) * I R 

DO 29 J - 1, IR 

L = J + LL 
M * J + KK 

29 DSUM = DSUM + CT(L)*CT(M) 

C ABOVE FORMAING T TRANSPOSE TIMES T. WHICH IS APPROX . OF B SHARP 
L = I + KK 

IF(QF) 31, 39, 32 

31 KK = KD3 + L 
D ( KK ) = D ( L ) 

EE ( L ) = DSUM 
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GO TO 30 

C-39 COMPUTE T TRANSPOSE * T ONLY.. PROVIDES INVERSE R SHARP 

39 EE( L ) = DSUM 

GO TO 30 

32 EE( L ) = D (L ) 

KK = KD3 + L 

0 ( KK ) = nSUM 
30 CONTINUE 

I F ( OF ) 41,66,41 

41 MV(1) = IR 

NV( 2 ) = IR 

C- LET P = 1ST MATRIX = EE,, 0 = 20 = D(K03+1) 

CCC ROLES OF P AND 0 ARE G I V F M TO 2 MATRICES AT 31, SWITCHED AT 32. 
COMPUTE D (K4) = 0*P, EE=D ( K 4 ) *0 , EE = EE — D { K 3 ) = 0*P*0 - 0 
C-FUNCTION IS NRM(0*P*0 -0)/NRM(Q) RESULT = A SCALAR 
CALL MULT ( 0 ( K03+1 ) , NV , FE , NV ,0 ( KD4+ 1 ) , MV ) 

NL = IR* IR 

CALL MULT ( 0 ( KD4 +1 ) , NV , P( KD3+ 1 ) , NV,EE,NV) 

DO 8 I = 1, ML 
KK = KD3 + I 

8 E E ( I ) = EE(I) - DC KK ) 

CALL NORM! FE ,NV,AM) 

CALL NORM ( D ( KD3+ 1 ) » NV , RR ) 

r.QUOTTFNT NFARS P.0 AS BSHRP APPROACHES THAT FITTING 2 MOOR-PENRSE AXIOM 
BDNRM = FN / FR 

9 RETURN 

66 BDNRM=FN 

C 66 IS A DUMMY REALLY WANT MAT RX MULT. ONLY. 

GO TO 9 

C SIDE COMPUTATIONS J W ANDREWS INF. SYSTEMS CO. MAY 1969 

END 





SUBROUTINE ANDRA (B , T , DPR , JP ) 

C- SUBROUTINE ANDRA DIAGONALIZES POS. DtF .SYMM. J ANDREWS I. S. CO. 

C - SUBR ANDRA CALLED BY PSEU J W ANDREWS, IMF. SYSTEMS CO. APRIL 1969 
IMPLICIT REAL*8 (D), INTEGERS (0) 

DOUBLE PRECISION B , T 

DIMENSION B I 400 ) , T(400), DPR ( 2 ) , JP ( 5 ) 

EOU I VALENCE I DP I ,FDI , IDD), (DCC, ICC ) , ( DMX , FMX ) , ( DR S , 1 1 S ) 

EQUIVALENCE ( DF Z 0, FZR 0, I Z ) 

DATA ICC, DFZO /Z40000000, O.DO/ 

C-DPR1 IS MAGNITUDE THAT IS CONSIDERED ZRO PIVOT MOST BE MO SMALER. 

C- DPR ( 2 ) IS TO RETURN FINAL PIVOT, SO THAT USER MAY TEST SMALLNESS. 
CC- ANDRA CAN BE USED ALL BY ITSELF TO GET INV., RANK OF POS SYMM. 

C NOT F THAT DSQRT HAS TO BE TAKEN OF PIVOTS ALONG THE DIAGONAL. 

C- MOTE I AM DELIBERATELY ALLOWING SOME PARAMETERS TO CHANGE ON SIJBSE- 
C-OUENT CALL DPR(l) CHANGES PIVOT SIZE A ROUGH TOLERANCE FOR ZRO. 

EF = SNGL I DPR 1 1 ) ) 

C- TEST- IS THIS AN I NT I AL I Z AT I ON CALL/ 

I F ( J P ( 2 ) ) 2, 1 , 2 
C INTIAL1ZE- FORM IDENTITY MATRIX 

1 OS = JP ( 4 ) 

ONT = OS*OS 

I F( OS .LT. 1 .OR. QNT .GT. 6400) GO TO 691 

DO 18 I = 1, ONT 
18 T ( I) = DFZO 

L_=_l 

0R1 = OS +1 
DO 1810 1 = 1, OS 
TIL) = 1 . DO 
1810 L = 0R1 + L 

DPR I 2 ) = DFZO 

C- SET RANK TO ZRO. TRIAL PIVOT VALUE TO ZERO. 

OKR = IZ 

C SET PIVOT CHOICE ITERATION AT 0 ALLOWANCE OF MO. ROWS+1 ITER. 
0 I TR = i Z 

2 CONTINUE 
200 CONTINUE 

C- 7 FR n OUT MAX DTAG. AND C T DTAG. TEMPORARY VARTARIFS 






FMX = FZRO 
I = IZ 

M £ IZ 

C- BELOW SEE IF ALL DIAG ELEMENTS TESTED YET 
L ■ 1 - OR1 

30 IF ( I «EQ. OS) GO TO 40 

1=1+1 
L = 0R1 + L 

DP I = B ( L ) 

C- GET CURRENT DIAG. ELEMENT FOR INTEGER, SINGLE PREC. TEST 
C- UPDATE L TO GET -NEXT- DIAG. ELEMENT 

C-BELOW TEST FDR DIG. ELEMENT DALREADYD REDUCED TP 1 . ( CODEMARKED) , ICC 
I F( IDD .EO. ICC) GO T 0 30 
I F ( FD I - FMX) 30,30, 32 

C- TEST FDR NEGLIGIBLE FL . PT. DTY. -TREAT THESE, AND MEG., AS ZEROS. 

32 IF C FD I .LT. EF) GO TO 30 

C- SET NEW MAX, 2BLE PREC., SAVE BEST ROW FOR PIVOT a OMRS) 

OMR = I 

DMX = DDI 
M = L 

GO TO 3 0 

40 CONTINUE 

400 I F ( M .EO. IZ) GD TO 505 " ' 

DRS = 1. DO /PSQRT(DMX) 

C SET INDEX OF FIRST ROW, OMR (PIVOT) COLUMN 

K = (OMR -1)*0S +1 

L = OMR 

DO 41 I = 1, OS 
DDM = b (L )*DRS 
T ( L ) = T(L)*DRS 
B ( L) = DDM 

C — SYMMETRICALLY, FORCE COLUMN TO SAME VALUE IN B ONLY 

B(K) = DDM 

K = K +1 

41 L = OS + L 

C FORCE PIVOT ELEMENT TO EXACT VALUE OF UNITY 

R(M) = 1.00 

C- NOW REDUCE ALL OTHER ROWS OF B, T, ELIMINATING COLUMN OF PIVOT VARIAB 
DO 460 I = 1, OS 
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O TEST FOR PIVOTAL ROW. OTHER ROWS 
I F ( I .EO. OMR ) GO TO 460 
COEFF, TO BE ZEROED CAM NOT RE PREVIOUS PIVUT. 

J = I - OS 
K = OS* I + J 

DRS = B ( K ) 

C RELOW TEST FOR A ROW ALREADY REDUCED, TO SKIP 

I F ( IIS .EO. ICC) GOTO 460 
C- GET COEFF IM PIVOT COLUMN TO BE ELIMINATED 
K = 0 M R * 0 S + J 
DMM = - B ( K ) 

L = OMR 

K = I 

DO 47 J =1, OS 

C- L IS ROW USED TO REDUCE, WITH PIVOT. 

C K IS CURRENT ROW THAT PIVOT GETS ELIMINATED FROM. 

B ( K ) = B ( K ) + B(L )*DMM 

T ( K) = T ( K ) + T(L)*DMM 

L = OS + L 
47 K = OS + K 

460 CONTINUE 

L = OMR 

DO 461 I =1 , OS 

C FORCE MOST OF PIVOT ROW TO ZERO. COMPLETES REDUCTION WITH 1 PIVOT/ 
B ( L ) = DFZO 

461 L = OS + L 

C FORCE PIVOT TO EiCODES) FOR ONE .. 

B ( M) = DCC 

C- SIGNAL NO LONGER FIRST TIME CALLED. 

JP( 2 ) = 1 

C — UPDATE EFFECTIVE RANK FOUND 
OKR = OKR +1 

DPR (2) = DM X 

J PI 5 ) = OMR 

C- NOW TEST -IS THIS AN ITERATION TO DO ONLY 1 ROW AT A TIME/ 

I F ( OKR . EO . OS) GO TO 480 

I F < JPC 1) .EO. IZ) GO TO 490 

C( AT THIS POINT, EITHER STOP WITH ONE ROW OR TRY NEXT. 

C HERE GET READY TO RETURN. RANK PARAMETER. 



480 JP( 3 ) = OKR 
RETURN 

C IF ENOUGH TRIES TO D ALL ROWS PLUS 1 MORE* QUIT. 

490 IF ( OITR .EO. 0R1) GO TO 480 
0 ITR = OITR +1 

GO TO 200 

691 CALL LNCNT(l) 

WRITE (6,1691) OS , QNT 

1691 FORMAT (' DIMENSION FRROR IN ANDRA NR = 1 » 1 4 « 5X . 1 NR#NC = 1 1 4 ) 
RETURN 

692 CALL LNCNT(l) 

WRITE (6,1692) 

1692 FORMAT (• ERROR IN ANDRA, FINDS NO PIVOTS') 

CHECK FOR DIAGONAL ALLOWING NO PIVOTS// 

505 I F ( J P ( 2 ) .EO. IZ .UR. OKR . GT, 0R1) GO TO 692 

GO TO 480 
END 


4 ^ 

-J 



oo 


SUBROUTINE DEC PM (A,B,C, 
C- SUBR OECOM FINDS 3 MATRICES 


B t C. ,E. JL 


dcm.kp. D) 

FROM WHICH USER 


CAN T,FT TTFrrWPTTTTTTnTT 


C INTO KROMECKER PRODUCT, POSSIBLY USING A SEPARATE PSEU CALL 

C.. SUBR. PECOM INFORMATION SYSTEMS CO. HAY 1969, J W ANDREWS 

IMPLICIT REAL-R (D), INTEGER*.? (0) 

DOUBLE P RFC I S I ON A,B,C,E, D, PIV 
COMMON /MAX/MAXRC 

DIMENSION A ( 400) ,B( ADO) ,C( 400) , E( 400) , 0(2000), 

2 JL (400) , DCM ( 3 ) , KP ( A ) , NV ( 2 ) , ND(2) 

EQUIVALENCE (NV, I), (PIV, ND( 1 ) ) , ( DFZO , 1 Z , 07, ) , (NP(2), J), 

1 ( N D ( 1 ) , OR 1 ) 

DATA OF Z 0 / 0 . DO / 

C- SET PARAMETERS TO CALL PSEU. TO FIND T TRANSFORMATION IN C. 

C — ASSUME INPUT A IS POS. DEFINITE SYMMETRIC 
OS = KP ( 3 ) 

QNT = 0 S * 0 S 

C-FRR EXIT IF ABSURD DIMENSIONS. 

I F ( OS .LT. 2 .OR. QNT .GT.MAXRC) GO TO 691 

C- SET COLUMN SIZE = RANK SIZE 

KP (4 ) = OS 
OL = K P ( 1 ) 

C- S ET SPECIAL PARAMS FOR PSEU CALL THESE ARE TO SUPPRESS THE WORK OF 
C RF-INVERTIMG PSFUPO INVERSE IN THE CASE WHERE A SINGULAR... 

K P ( 1 ) = - KP ( 1) -1 

C- CALL PSEU P TO GET MATRIX T. IN C 

CNOTE THE LAST 3 MATRICES OF THE 5 IN D USED ONLY IF PSEUP 0ITERATES5) 
I = KP( 2) 

CALL PSEU P( A.B.C.E.DCM.KP, D) 

KP< 1 ) = OL 

IF ( I .ME . IZ ) GO TO 38 

C/ PLEASE DO NOT TRY TO TAKE A.S.P. NAMES FOR MATRICES HERE. 

C- SUCH MATRICES WERE NOT RETURNED BY ASP, NOR BY MY. ROUTINE. 

C 

on 13 I = 1, QNT 

13 D ( I ) = C ( I ) 

NV ( 1 )=0S 

NV ( 2 ) =0S 

N D ( 1 ) = 2 

C- MD IS PART OF FLAG (PIV) RETURNED BY INV. 

CALL INV(D,NV,PIV, JL) 



— fvni— 


C- TEST T D SEE IF PIV IS ZERO = ERR, MATRIX NOT INVERTIBLE. 

I F ( N 0 ( 1 ) .EO. IZ ) GO TO 692 

38 CONTINUE 

KU - ONT 

C- P IS TO HOLD PERMUTATION MATRIX SUCH THAT THAT P*BE* PTR = ER 

C- CR has ALL ONES MOVED TO EXTREME UPPER LEFT OF DIAGONAL. 

0* NOW SET UP TO MAKE P PERMUTATION MATRIX P = D(KO +1) 

OR 1 = OS +1 

C- ZERO OUT P, WILL BE ZEROS AND ONES 

DO 39 I = 1, ONT 

C-ZERO HOUSEKEEPING ARRAY ONLY NEED FIRST COLUMN. 

JL( I ) = IZ 
K = KD + I 

39 D ( K) = DFZO 
L =1 

r/j =j 

OL =1 

DO 780 K = 1, OS 

I F ( R( L) ) 7803, 7801, 7803 
7803 J = JL(OL) 

CHECK FOR ROM OF DIAG. THAT NEEDS A 1 MOVED INTO IT 

I F ( J .EO. IZ) GO TO 78fe 
I = (K-1)*0S + J + KD 

C-PIJT 1 IN P TO MOVE K.K TO J.K POSITION (2 PFRMUTAT IONS ) P,, P TRANS 
C* THE EFFECT IS TO MOVE K , K TO J,K/ THENCE TO J,J. 

D ( I) = 1 . DO 

O L = OL +1 

C/MARK THIS 1 AS ZRO TO BE FILLED— IT IS MOVED IjP AND OIJTOF HERE 

7801 JL ( M ) = K 

M = M +1 

GO TO 780 

786 J = KD + L 

C .MAKE PART OF IDENTITY AT 786 DON-T NEED TO MOVE 1 TO A HOLE. 

D(J) = 1 . DO 
780 L = 0R1 + L 

C RETURN. MATRICES COMPLETED E WITH IR IAS DELIBERATELY LEFT OUT. 

RETURN 

691 CALL LNCNT(l) 

WRITE (6.1691) OS. ONT 
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1691 FORMAT (• DIMENSION ERROR IN DECOM NC= 1 , 1 A , 5 X , ' NR* NC = ' » I A ) 
RETURN 

692 CALL LNCNT(l) 

WRITE (6,1692) 

1692 FORMAT (• ERROR IN DECOM PIVOT=ZERO') 

KP( 4)=-QS 

GO TO 38 
END 




APPENDIX C 


USE OF VASP ON AMES’ TSS 
NONCONVERSATIONAL (BATCH OR RJE) 


In using VASP on TSS, the system must be told about the job library in which the VASP 
subroutines are located, the source of input data, and the location to send output data; and the 
block data program must be loaded. 

A procedure has been written for doing this automatically. The call to the procedure is 

VASP$$ [input data set] [, output data set] 

The procedure will then perform the steps indicated above. If the first parameter is omitted, the 
data will be taken from SYSIN, which is from cards in your data deck. If an input data set is 
named, then the data will be taken from the named data set, which must have been stored 
previously. 

Likewise, if the second parameter is omitted, the output will be placed in SYSOUT, for 
printing on the high-speed printer. If an output data set is named, the output will be placed in 
that data set. 

If the name of the input or output data set must be changed, use the procedure call 

CHNGIN [new input data set name] 

CHNGOUT [new output data set name] 

These two procedures will then change the DDEF to the new data set name. If the 
parameter is omitted, the new data set name will be SYSIN or SYSOUT. A listing of these 
procedures is included in this appendix. 


CONVERSATIONAL 


Provisions have also been made to allow conversational use of the VASP program, so that the 
user can easily perform matrix operations. The operations can be strung together in a sequence as 
desired with as much output as desired. The user indicates the operations by use of Fortran 
statements, and may not only call the VASP subroutines, but also may execute any other Fortran 
statements that he wishes. 

Data are requested for the program by means of subroutine INPUT, allowing free-form data 
from the typewriter. If Fortran type input is used, the data should also be obtained from the 
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typewriter. If you try to use an input data set, INPUT will also read the same data set. 
Variables may also be set by Fortran arithmetic statements. 

Output may be from the VASP subroutine PRNT, or any Fortran WRITE or PRINT 
statement. Two standard formats are available if desired for unlabeled output. 

The program automatically dimensions 14 arrays to the desired size, and the user may 
supply his own names to 7 of them. 


Usage 

The use of conversational VASP is demonstrated by the accompanying figure (fig. 9). 

Lower case letters are input and upper case are the computer responses. Detailed comments on 
the various statements follow. To start, the user calls VASP$$ (line 1) as for nonconversational 
usage. If desired, an output data set may be named. Line 2 lists the DDNAME being used. 

The next two lines (lines 3 and 4) indicate where input and output are to reside. The 
computer then gives an underscore, after which the procedure “CONVASP” is called. The param- 
eters of this procedure are first the total number of elements in a matrix, followed by up to 
seven matrix names. If the parameters are defaulted, the system will select matrices with 9 ele- 
ments, and name the matrices A, B, C, W, X, Y, Z. In addition, 7 dummy matrices D1 through 
D7 are available for use. In the figure, all matrices are to be dimensioned 16 (line 5), the second 
matrix is to be renamed F, and the Z matrix is to be renamed FSTAR. That is, if you wish 
to rename a specific matrix, put a dollar sign in front of the original name and then equate it to 
the desired name as in the example. Fourteen arrays, NA through ND7, used for dimension 
information, are also defined and renamed to agree with the working matrices. 

Lines 6, 7, and 8 then define the matrices available. Note that no 1- element variables are 
defined. The user may define them in his program but they will not be available from one 
computation to the next. 

The computer will then ask for FORTRAN STATEMENTS?. At this point, a data set 
SOURCE. MNPGSS has been set up for editing and the necessary DIMENSION and other initial- 
izing statements have been stored. These statements are listed in figure 12, lines 4600 through 
6000. The computer prompts the user with 100 and the user may enter any Fortran statements 
he wishes. The full power of the text editor is available at this point. 

In the example, we have entered four statements, lines 10 through 13. Note that we have 
defined a single variable t for use in the etphi statement. The value of this variable will not 
be remembered by the system. 

After completing the desired Fortran statements, the user requests compilation by entering 
_CMPL (line 14). The computer then indicates that compilation is proceeding (line 15) and will 
give the usual error messages if the compile is unsuccessful. After compilation the program is 
automatically executed, and the first item in the execution is a request for data from the INPUT 
subroutine (line 16). Data are entered free style as in line 17, with the elements of the matrices 
152 
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1 

_ 2 . 

3 

4 

5 


vasp$$ 

DDNAME=JBLB0001 
INPUT FROM TERMINAL 
OUTPUT TO TERMINAL 
£onvasp 16,,f,$z=fstar 


“9 

10 

11 


TZ 

13 


8 *****MATftlCES AVAILABLE, ALL DIMENSIONED 16, ARE; 

7 A,F,C,W,X,Y,FSTAR; FOR INPUT OR COMPUTATIONS 

8 01,02,03,04,05,06,07; FOR COMPUTATIONS O NLY 
*****TORTRAN STATEMENTS? 

0000100 t=l. 0 

0000200 call etphi (f ,nf . t,a,na,dl, 32) 

(f , nf , f _r ;i) 


0000300" 

0000400 


cal 1 prnt 
cal 1 prnt (a, na, 'a 


M> 


-T5~ 

*****MNPG$ip NOW COMPILING***** 




16 

DATA? 




17 

f-1. 1,1. 2, 1.3, 1.4, 2. 1,2. 2, 2. 3, 2. 

4, 3. 1.3. 2. 3. 3. 3. 4 

,4. 1,4. 2,5.3 

, 4. 4,nf =4, 4* 

18 

F MATRIX 4 ROWS 

4 COLUMNS 



19 

1. 1000000D 00 2 . 1000000D 00 

3. 1000000D QO 

4 . 1 0000QQD- 

DO— - 

" 2 ir 

"T720(T0000D 00 2 . 2000000D 00 

3. 2000000D 00 

4. 2000000D 

00 

21 

1 . 3000000D 00 2. 3000000D 00 

3 . 3000000D 00 

4. 3000000D 

00 

22 

1. 4000000D 00 2.4000000D 00 

3 . 4000000D 00 

u. tifinnnnnn 

DO 

23 

A MATRIX 4 ROWS 

4 COLUMNS 



24 

7. 72 5164 7D 03 1.3690166D 04 

_ 1^9656.1670 04_ ... 

7 .JI 6221680 04 - .... 

75 

8.0162773D 03 1.4208825D 04 

2. 0399372D 04 

2. 6590919D 

04 

26 

8. 3083898D 03 1.4725483D 04 

2. 1143577D 04 

2. 7559671D 

04 

27 

8. 6005023D 03 1.5243142D 04 

2. 1885782D 04 

2.85294220 

D4 

28“ 

*****Cf)MPUTI NG DONE***** 




29 

recmpt 




30 

DATA? 




~3T 

f =1.11, "1 722717 33 




32 

f(4)=2. 11,2. 22,2. 33, 




33 

f(7)-3. 11,3. 22,3. 33 




“34“ 

~"nt =3, 3 




35 

* 





Figure 9.— Example of conversational VASP. 




36 

F MATRIX 

3 ROWS 

3 COLUMNS 

37 

1. 1100000D 

00 

2. 1100000D 

00 

3.1100000D 

00 

38 

1.2200000D 

00 

2. 2200000D 

00 

3u22.CLQQQ.OD 

..OIL 

39 - 

" ' 1. 3300O(TOD 

00 

2 . 3300000D 

00 

3.3300000D 

00 

40 

A MATRIX 

3 ROWS _ 

3 COLUMNS- 

tIT- 

l.'5033"4T5D 02 

2 . 686661 ID 

02 

3. 8799809D 

02 

42 

1.57051530 

02 

2. 8346117D 

02 

4. 0787080D 

02 

43 

1. 6476893D 

02 

2. 9625623D 

02 

4. 2874352D 


"44 

***COMPUT 1 NG" 

TJUTTE**-* 





45 Jiewrt _ . 

"46 ***FORTRAN STATEMENTS? 

47 0000100 call mult (a,r>a # x,nx/ y/ny) 

48 0000200 call prnt (y,'y ',1) 

"49 OUUU JU0_rev i se 2W 

50 0000200 call prnt (y/ny,'y ’,1) 

51 0000300_1 nser t 150 _ 

52 0 O' 0TJY5 0 p FT rt"~6 7f 

53 £mpl 

54 *****MNPG$$ NOW COMPILING***** 

“55r-D7n7T? 

56 nx=3, 1, x=l. 0, 0. , 0. * 

57 1 . 1 10000D 00 1. 22000 0P 0 0 1.330000D_00 2.1100000 00 2. 2200000 00 2 . 33nnimD-.Ql) 

-58- 3 . 11D OMtT 0 tr 3.2 2 0 0 0 OD 00 3.3300066 66 3. 200000D 00 3.300000D 00 3.400000D 00 

59 4. 100000D 00 4. 2000000 00 4.300000D 00 4.400000D 00 

' 6e 7 MATRIX T~ ROWS 1 COLUMNS 

61 1 . 5033413D 02 

62 1. 5705153D 02 __ 

63 I76476T33D 02 

64 *****COMPUTI NG DONE***** 

65 _rewr t 150 — 

66 urruimro catl mult (a,na,x,nx,y,ny) 

67 0000150 call add (x, nx, y, ny, w, nw) 

68 0Q00250 call p rnt (w.nw/w Vl) 

"69 000 035 0 call prnt x,nx, 1 x 1 , 1 ) 

70 0000450_cmpl 

71 *****MNPG$$ NOW COMPILING***** — — 

Figure 9.- Example of conversational VASP - Continued. 



IT 

73 

74 

0000350 T **"* r NOT FOUND WHERE 
0000350 CALL PRNT X/NX/'X ',1) 
#350/ call prnt (x,nx,'x M) 

REQUIRED 



75" 

76 

77 

MODIFICATIONS? 

n 




78" 

79 

DATA? 

* 




-80" 

81 

82 

W MATRIX 3 ROWS 

1. 5133413D 02 
1. 5705153D 02 

1 COLUMNS 



■83" 

1.6476893D 02 




84 

X MATRIX 3 ROWS 

1 COLUMNS 



85“ 

86 

87 

~T7inr<JM0OD 00 
0.0 
0.0 




■ir 

89 

90 

*****COMPUTING DONE****'*" ~ 

logoff 

9.729 CPU SECONDS/ 05/11/71 AT 

11:35/ FST04 

R3984 



Figure 9.— Example of conversational VASP — Concluded. 
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being entered columnwise. Do not forget to input the matrix dimensions such as NF in the 
example. Data entry is ended with an *. Execution of the program continues; lines 18 through 
27 display the requested output, and line 28 indicates completion. 

At this point (line 29) the computer gives an underscore and the user may do anything he 
wishes. In the example, we are going to recompute with the same program, using new data. 
Accordingly, the user asks for RECMPT (line 29). The program is again executed, and new data 
are asked for (line 30). They are entered in lines 31 through 35, using a different style than in 
line 17 to show the flexibility available. On completion of the data entry, the results are 
printed in lines 36 through 44. 

At this point, it is desired to rewrite the entire program, so the user issues the command 
REWRT (line 45). The system, as at line 9, prompts the user with “FORTRAN STATEMENTS,” 
and a line number (lines 46 and 47), after which the user enters Fortran statements as desired. 

In >the example, line 48 is entered incorrectly and then corrected (lines 49 and 50). Following 
this, a line 150 was inserted (lines 51 and 52). Then a CMPL was issued (line 53) to compile 
and execute the program. New data were entered at line 56, and lines 57 through 59 are the 
output requested by the statement “print 6,f.” Note that all 16 elements of f are printed 
' using one of the two FORMAT statements compiled into the program for convenience (see 
lines 5900 and 6000 of VASPPROC, fig. 12): 

6 FORMAT (IX, 1P6D 13.6) 

13 FORMAT (IX, 1P4D20. 13) 

These statements request the output of a 6 decimal number or a 1 3 decimal number. In the 
example, we are printing a 6 decimal number. The remainder of the output is then printed 
(lines 60 through 63). 

Now, it is desired to rewrite only a portion of the program from line 150 on. Accordingly, 
the REWRT command is issued with a parameter (line 65). The system then erases 
SOURCE.MNPGSS from line 150 inclusive to the end. It then lists that portion of the program 
being used, in this case, line 100 only (line 66) and prompts the user for additional lines with a 
line number (line 67). The user then adds lines as desired (lines 67 through 69) and requests a 
compile (line 70). It can be seen that line 69 is missing a left parenthesis so the compiler prints 
a diagnostic and requests the line be corrected (lines 72 and 73). The correction is entered 
(line 74), after which the compilation is completed (lines 75 through 77). No data are needed, 
so the data request (line 78) is answered with * only (line 79). The results are printed on lines 
80 through 88. Since no more computations were desired, a Logoff command was issued 
(lines 89 and 90). 


Housecleaning 

A procedure called “CLRVASP” is available. This procedure erases all data sets that have 
been set up by the various other procedures, and allows the user to keep his storage low. Use of 
the routine is not required since the other procedures have appropriate erase statements as 
needed. 
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LISTINGS AND FLOWCHARTS 


Figure 1 0 shows all the procedures associated with VASP, and indicates what each one does. 
A complete listing of the procedures is given in figure 1 1 . Figure 1 2 is a listing of data set 
VASPPROC. If the user executes this data set, it will generate all the procedures and place them 
in the user’s USERLIB. 


TSS ACCESS 


For access to the VASP program, an Ames TSS user should issue the following statements: 

SHARE VASP, FSTJSW, VASP 

which allows access to the VASP subroutines 

SHARE VASPPROC, FSTJSW, VASPPROC 
EXECUTE VASPPROC 

which first allows access to a data set containing the various procedures, and then enters these 
procedures in the user’s USERLIB. Note that the EXECUTE command sets up a batch job, and 
that the procedures will not be available until that batch job is completed, and the user has 
issued either a LOGOFF or ABEND command. After once issuing these commands the user 
need only call the procedure, as discussed earlier. 

Further, for conversational use, issue the command 

SHARE VASP1, FSTJSW, VASP1 

which allows access to the proper version of subroutine INPUT. 
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VASP$$ 

JBLB VASP 
LOAD BLKDTA$$ 
Input & Output DDEF 
Default Options 


CONVASP 
JBLB VASP1 
DISPLAY Matrix Names 
Edit SOURCE. VASPMN$$ 

Compile VASPMNSS 
Load VASPMN$$ 

Edit SOURCE. MNPGSS 
EXCERPT beginning of Fortran 
Programs 

Display FORTRAN STATEMENTS? 


CHNGIN 

Change Input DDEF 


CHNGOUT 
Change Output DDEF 


CMPL 

Add end of Fortran Program 
Display MNPGSS Now Compiling 
Compile MNPGSS 
Call MNPGSS 

Display COMPUTING DONE 


CL R VASP 
Erase all Programs & 

Data used by Conversational 
VASP 


RECMPT 
Call MNPGSS 


REWRT 

EDIT SOURCE. MNPGSS 
EXCISE Program 

Display FORTRAN STATEMENTS 


REWRT N 

Edit SOURCE. MNPGSS 
EXCISE from Statement N 
to last 

List program 


Figure 10.— Flowchart VASP procedures. 
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CHNGIN' 0000000 PROCDEF CHNGIN ~ 
CHNGIN 0000100 PARAM $ INPUT 
CHNGIN 0000200RELEASE FT05F001 
ettNGIft— OW1WW— 1-F-^ ♦ttNPUT*- ; D 0fTF~F 
CHNGIN 0000400 IF '$INPUT' = ,, ;DISPLAY 


CHNGOUT 0000000 PROCDEF CHNGOUT 
CHNGOUT 0000100 PARAM {OUTPUT 

-chngolhf -ewfim- pehase- ftoc-fooi - - 

CHNGOUT 0000300 IF '{OUTPUT' ' ;DDEF FT06F001, , {OUTPUT; D I SPLAY 'OUTPUT PLACED IN DATA SET {OUTPUT 1 
CHNGOUT 0000400 IF '{OUTPUT* =";DI SPLAY 'OUTPUT TO TERMINAL' 


0000900 

^RYASPMOoocrtri^cnprTn.RVASP 

CLRVASP 0000050 END 
CLRVASP 0000100 UNLOAD MNPG{{ 

CLRVASP ”0000200 UNLOAD VASPMN${ 

CLRVASP 0000300 ERASE SOURCE . VASPMN{{, SOURCE. MNPG{{, USERL I B (VASPMN{{ ) 

CLRVASP 0000400 ERASE USERL I B (MNPG{{ ) 

-etRVAST> 00005150 RELEASE"' VAST! 

CLRVASP 0000600 DISPLAY ’*****ALL CONVERSATIONAL VASP PROGRAMS CLEARED***** 1 
CLRVASP 0000700 DISPLAY '*****YOU MAY RESTART WITH COMVASP******** ' 


CMPL 0000000 PROCDEF CMPL 

UMPL 0000030 D EF AULT - SYS I NX =E " — — 

CMPL 0000050 EDIT SOURCE. MNPG{{ 

CMPL 0000100_EXCERPT SOURCE. VASPMN{{, , 1600, 1700 

CMPL 0000150 END 

CMPL 0000200 DISPLAY ' *****MNPG{$ NOW COMPILING*****' 

CMPL 0000220 DEFAULT LIMEN=N 

CMPL 0000250 UNLOAD MNPG{{ “ 

CMPL 0000270 ERASE USERL I B (MNPG{{ ) 

CMPL 0000300 FTN MNPG{{,Y 

CMPL 0000400 LOAD MNPG{{ 

CMPL 0000600 CALL MNPG${ 

CMPL 0000700 DEFAULT LIMEN=W 

CMPL 0000800 01 SPLAY f *****COMPUTING DONE*****' 


Figure 11.— List of VASP procedures. 
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X0NVASP D0OOO0O PROCDEF CONVASP 
CONVASP 0000020 PARAM $N,$A,$B,$C,$W,$X,$Y,$Z 

CONVASP 0000040 DDEF VASP1, VP, VASPl / OPTION=- l OBLI B 

- eoNVA S P ooo ooso -job-libs -sysultb - 


CONVASP 0000110 
CONVASP 0000140 


C O N V AS ~p— 000 OT7 0 01 SPLAT 

CONVASP 0000200 DEFAULT SYSINX=E 
CONVASP 0000250 DEFAULT LIMEN=N 

— c onvas p - oooosoo— Eurr source; vas p m ntt 

CONVASP 0000340_EXC I SE 1, LAST 
CONVASP 0000380 INSERT 100 

CONVASP" 0000400” IMPLICIT REA 

CONVASP 0000430 COMMON /ASP/ 

CONVASP 0000460 


DISPLAY ' *****MATR I CES AVAILABLE, ALL DIMENSIONED $N, ARE;' 
DISPLAY ' $A,$B,$C,$W,$X,$Y,$Z; FOR INPUT OR COMPUTATI 

"01 SPLAT Dr / U2;UJ / D4/D3 / D6 / DT;' FOTT COMPUTATIONS ONLY' 


CONVASP 


IMPLICIT REAL*8CA-H”,0-ZT ‘ 

COMMON /ASP/ $A($N),$B($N),$C($N),$W($N),$X($N),- 

1 $Y($N),$Z($N),D1($N),D2($N),D3($N),D4($N),D5($N) 

- 2 - n$a( 7 ) , ct2 1 , N$wm, N$XTTT,- 

3 N$Y(2),N$Z(2),ND1(2),ND2(2),ND3(2),ND4(2),ND5(2) 
COMMON /MAX/ MAXRC 

MAXRC=$N 

10 PRINT 15 
15 FORMAT C DATA?') 

1 '$X' < $X/$Y' / $Y,'$Z',$Z,'N$A' / N$A # - 

2 , N$B' / N$B / , N$C' / N$C / , N$W',N$W, 'N$X',N$X,-_ 

- -• 3 , N$Y , ,N$Y/'N$Z , / N$'Z)“ 

13 FORMAT ( IX, 1P4D20. 13 ) 

6 FORMAT (IX, 1P6D13 . 6 ) 


CONVASP 0000520 3 N$Y(2),N$ 

CONVASP 0000550 COMMON /MAX/ MAXRC 
-COWVAS-P-OOOO500 -MAXRC=$N 
CONVASP 0000700 10 PRINT 15 

CONVASP 0000800 15 FORMAT C DATA?') 

-CONVASR 0 0 00 84 0 — CAT t~T NPUT ( A’ , $ A', 
CONVASP 0000880 1 '$X',$X,'$Y 

CONVASP 0000920 2 'N$B',N$B, ' 

-convasp oooooso - -■ 3 , i»$y , ,n$y/* 

CONVASP 0001000 13 FORMAT ( IX, 1P4D20. 13 

CONVASP 0001100 6 FORMAT ( IX, 1P6D13 . 6 ) 

ctJNVA^p- -ooottoo — return — 

CONVASP 0001300 END 
CONVASP 0001400_END 

CONVASP 0001450 ERASE USERL IBCVAS PMN$$ ) 
CONVASP 0001500 FTN VASPMN$$,Y 
CONVASP 0001600 XLIST VASPMN$$ 


, D6($N),D7($N) , ; 
, ND 6(2), ND 7(2) 


CONVASP 0001600 
— CONVA SP" 00 0170 tr 
CONVASP 0001770 


DEFAULT LI MEN=N 
EDIT SOURCE. MNPfi$$ 


CONVASP 0001800 EDIT SOURCE. MNPfi$$ 

"CONVASP 0OO10OO_TXCTSri,LAST 
CONVASP 0001950 INSERT 1,1 

CONVASP 0002000_EXCERPT SOURCE. VASPMN$$, , 100, 1500 
-CONVA SP 00 0 21 00 DTS PLAY - • *****FTJRTRAN' STATEMENTS^' 
CONVASP 0002150 DEFAULT LIMEN =W 
CONVASP 0002200 INSERT 100 


Figure 11.— List of VASP procedures — Continued. 



BECMPT 0000000 PROCDEF RECMPT 
RECMPT 0000100 DEFAULT LIMEN =N 
RECMPT 0000200 CALL MNPG$$ 

RECMPT 0000300 DEFAULT LIMEN-W 

RECMPT 0000400 DISPLAY 1 ***COMPUTI NG DONE***' 


REWRT 

REWRT 

REWRT 

REWRT 

REWRT 

REWRT 

REWRT 

REWRT 

REWRT 

REWRT 


0000000 PROCDEF REWRT 
0000100 PARAM $ L I NE 
0000200 DEFAULT LIMEN -W 
0000400 DEFAULT SYSINX=E 
0000500 EDIT SOURCE. MNPG$$ 

G0006G0_EXCISE $ L I NE, LAST 

0000700 IF ' $LI NE ' 100 ' ;D I SPLAY 1 ***FORTRAN STATEMENTS?' 
0000800 IF '$LINE ,-, = '100';LIST 100,LAST 
0000900 DEFAULT SYSINX-G 
0001000 INSERT $LI NE 


^ASP$$ 

VASP$$ 

VASP$$ 

VASP$$ 

-VASPU- 

VASP$$ 

VASP$$ 

V ASP £$ 
VASP$$ 


0000000 PROCDEF VASP$$ 

0000100 PARAM $ I NPUT, $OUTPUT 

0000150 DEFAULT $N=9,$A=A,$B=B / $C=C,$W=W / $X=X / $Y=Y,$Z=Z,$LI NE=100 
0000200 JBLB VASP 
0000300LOAD BLKDTA£$ 

0000400 IF ' $ I NPUT' ^=";DDEF FT05F001, , $ I NPUT;D I SPLAY 'INPUT FROM DATA SET $INPUT' 

0000500 IF ' $ I NPUT ' =";DISPLAY 'INPUT FROM TERMINAL' 

0000600 IF ' $OUTPUT ' " , -";DDEF FT06F001,,$OUTPUT;DISPLAY 'OUTPUT PLACED IN DATA SET $OUTPUT' 
0000700 IF ' $OUTPUT' =";DISPLAY 'OUTPUT TO TERMINAL' 


Figure 11.— List of VASP procedures - Concluded. 



~ 00030 LOGON USERID, ,9 

00060 PROCDEF CHNG IN - - ' ' ~ 

00090_EXC I SE 1, LAST 
00100 PROCDEF CHNGIN 

00200— P A RA M - -$ f NPUO 

00300RELEASE FT05F001 

00400 IF 1 $ I NPUT ' "’= ' ' ;DDEF FT05F001,,$INPUT;DI SPLAY 'INPUT FROM DATA SET $INPUT' 

00500 IF '$1NPUT L ;DLSPLAY 'INPUT FROM TERMINAL* 

00540_ PROCDEF CHNGOUT 
005 80_EXC I SE 1, LAST 

00600- PROOBEF-CHNGeUT- - ■ - - 

00700 PARAM $OUTPUT 
00800 RELEASE FT06F001 

009-00- - I F 1 $-OOTPUT l * 1 = 1, ;DDEF- FT06F001, , $GUTPUT;D+SPtAY 'OUTPUT ^ACEtT 1 N DATA OET ^OUTPUT* 
01000 IF 1 $OUTPUT' =";DI SPLAY 'OUTPUT TO TERMINAL' 

0 1 0 4 0_ PROCDEF CLRVASP 

OLOOO^EXOfS-E — 3r>-LA5T 

01100 PROCDEF CLRVASP 
01200 END 

01300 UNLOAD MttPGtt - - — - — ' 

01400 UNLOAD VASPMN$$ 

01500 ERASE SOURCE. VASPMN$$, SOURCE .MNPG$$, USERL I B ( VAS PMN$$ ) 

00600 ERASE US€fttt-B ( NTO P Ot»-) 

01700 RELEASE VASP1 

01800 DISPLAY '*****ALL CONVERSATIONAL VASP PROGRAMS CLEARED*****' 

01900 D ISPLAY '*****YOU MAY RESTART WITH CONVASP*™*****' 

0 1 9 4 0_ PROCDEF CMPL 
01980_EXC I SE 1, LAST 

02000 PR OC D EF - CMPL— 

02100 DEFAULT SYSINX=E 
02200 EDIT SOURCE. MNPG$$ 

02300 EXCERPT SOURCE. VASPMN$$, , 1600, 1700 - - ~ 

02400 END 

02500 DISPLAY '*****MNPG$$ NOW COMPILING*****' 

02600- BE-FAtttT- LLMEtDTt 

02700 UNLOAD MNPG$$ 

02800 ERASE USERL I B (MNPG$$ ) 

02900 FTN MNPG$$,Y ' 

03000 LOAD MNPG$$ 


Figure 12.- List of data set VASPPROC. 



03100 CALL MNPG$$ 

•9 -32 - fl 0 — DEFAULT L I MEN g W 

03300 DISPLAY » *****COMPUTI NG DONE***** 1 
0334 0 PROCDEF CONVASP 

0 3 3fttHE X€+S E - ±, t AST 5 

03400 PROCDEF CONVASP 

03500 PARAM $N, $A, $B, $C, $W, $X, $ Y, $Z 

03600 — B DEF VASP1, VP, VA& P1 > OPT! O N -JO B L LB 

03700 JOBLIBS SYSULIB 

33*60- - -D-f^^AY— f*-»***M A TfH-eES- AVA-HABLTi / — A Lt-D -m ENStO UE-B-tfL , A -R E; - 1 - 

03900 DISPLAY $A, $B, $ C/ $w, $ X, $ Y, $Z; FOR INPUT OR COMPUTATIONS' 

04000 DISPLAY ' 01,02,03,04,05,06,07; FOR COMPUTATIONS ONLY ' 

0 4 100 — DEF A ULT SY 54 N - X ~ E 

04200 DEFAULT LIMEN=N 
04300 EDIT SOURCE. VASPMN$$ 

frWQ ^CXC -t- SE 1, t A ST — — 

04500 INSERT 100 

04600 IMPLICIT REAL*8 (A-H,0-Z ) 

WOO COMMO N /ASP/ — $A($N) , $ D ( $ N3 , $C( $N F ,$W( $N) , $X( $N) , 

04800 1 $Y($N),$Z($N),D1($N),D2($N),D3($N),D4($N),D5($N),D6($N),D7($N),- 

04900 2 N$A(2),N$B(2),N$C(2),N$W(2),N$X(2),- 

3-tttYffr}-, N$Z( -2-^ , ND Iffl, ND2 ( 2) , ND3 ( 2 ), H D4 ( 2 ) , ND5 ( 2 ) , N D 6 ( 2 ) , ND 7tT? — 

05100 COMMON /MAX/ MAXRC 
05200 MAXRC=$N 

03 3 6 0 — 10 PR I NT 15 

05400 15 FORMAT (' DATA?') 

05500 CALL INPUT ( ' $A' , $A, ' $B ' ,$B, ' $C ' , $C, ' $W' , $W, - 

05600 1 ~ 1 $X * , ~$X ^ ' $ Y ' , $ ,tt7 4 ~ N$A l , N$A, ^ 

05700 2 ' N$B ' , N$B, ' N$C ' , N$C, ' N$W' , N$W, ' N$X ' , N$X, - 

05800 3 'N$Y',N$Y, 'N$Z',N$Z) 

05900 — 13 FORMAT (IX, 1 P 4 P 20. 13) 

06000 6 FORMAT ( IX, 1P6D13. 6) 

06100 RETURN 

6*200 ~ENfr 

06300 END 

06400 ERASE USERLI B( VASPMN$$ ) 

0 65- 00 — FF N VAS PM N $ $, V 

06600 XL I ST VASPMN$$ 

06700 LOAD VASPMN$$ 


Figure 12.- List of data set VASPP ROC -Continued. 
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OfifrDO— DEFAULT t+MEN=W — — 

06900 EDIT SOURCE. MNPG$$ 

07000 EXCISE 1,LAST 

07100 — I NSERT 1,1 

07200 EXCERPT SOURCE . VASPMN$$, / 100 / 1500 

07300 DISPLAY ' *****FORTRAN STATEMENTS?' 

DEFAutT-tmeN-^w 

07500 INSERT 100 
0 7 5 4 0_ PROCDEF RECMPT 

075 8 0_EXC I SE" f ; LAS T 

07600 PROCDEF RECMPT 
07700 DEFAULT LIMEN »N 

8 7800 ~CAlt M NPG$$ “ — 

07900 DEFAULT L I MEN=W 

08000 DISPLAY ' ***COMPUTI NG DONE***' 

0- 8 -0-4 0_ - PR0CD CF REWRT 

08080_EXC I SE 1, LAST 
08100 PROCDEF REWRT 
08200 PARAM $ L I NE 

D 8 30 6 — DE FA U LT L IM EN =W — 

08400 DEFAULT SYSINX=E 
08500 EDIT SOURCE. MNPG$$ 

HQ CA fi C V A I r C ± li r — I A C T - - — — — 

U obUU EXCISE $ Lt NE/ LAST 

08700 IF *$LI NE* = * 100* ;D1 SPLAY ' ***FORTRAN STATEMENTS?' 

08800 IF '$LINE'“ , = '100';LIST 100,LAST 

08 9 00 — DEFAULT S ' YStffl fr - e 

09000 INSERT $LI NE 
09040_ PROCDEF VASP$$ 

6 9 8fr fl_EXC I SE 1, LA S'? 

09100 PROCDEF VASP$$ 

09200 PARAM $ I NPUT, $OUTPUT 

09300 — DEFAULT $N=9, $A=A, $B = B, $ C=C , $W=W, $ X= X , $ Y=Y, $ Z°Z, $ L I N E = 100 

09400 JBLB VASP 
09500LOAD BLKDTA$$ 

08800 TF“ 1 $ I NPUT' -^= ++ 7tJnFFFT05FOOl / , $ I NPUT; D I SPLAY ' TNPUT^FROM - D" AT A" SET $ I NPUT' 

09700 IF ' $ I NPUT' =";DISPLAY 'INPUT FROM TERMINAL' 

09800 IF ' $OUTPUT' " ;DDEF FT06F001, , $OUTPUT;D I SPLAY 'OUTPUT PLACED IN DATA SET $OUTPUT' 
09900 IF 1 $O t J TP UT' ' g " ; DI SPLA Y 'O U TP U T TO TER M I NAL' 


Figure 12.- List of data set VASPPROC - Continued. 



1000 0 L I ST 

10100 END 

10200 — TALK * 

10300 VASP PROCEDURES NOW READY. DO 'ABEND' TO MAKE THEM AVAILABLE 
10400 LOGOFF 


Figure 12.- List of data set VASPPROC - Concluded. 
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